Scalable arrays of micro-Penning traps for quantum computing and simulation

TL;DR

This paper proposes a scalable 2D array of micro-Penning traps for quantum computing and simulation, enabling dense ion arrangements with high coupling rates and local control, overcoming power limitations of RF traps.

Contribution

It introduces a novel scalable architecture of micro-Penning trap arrays with stable ion motion and demonstrates feasibility for quantum simulation and computing.

Findings

Dense ion arrays with high coupling rates are achievable.

Stable ion motion is maintained despite trap imperfections.

Local control enables high-fidelity quantum gates.

Abstract

We propose the use of 2-dimensional Penning trap arrays as a scalable platform for quantum simulation and quantum computing with trapped atomic ions. This approach involves placing arrays of micro-structured electrodes defining static electric quadrupole sites in a magnetic field, with single ions trapped at each site and coupled to neighbors via the Coulomb interaction. We solve for the normal modes of ion motion in such arrays, and derive a generalized multi-ion invariance theorem for stable motion even in the presence of trap imperfections. We use these techniques to investigate the feasibility of quantum simulation and quantum computation in fixed ion lattices. In homogeneous arrays, we show that sufficiently dense arrays are achievable, with axial, magnetron and cyclotron motions exhibiting inter-ion dipolar coupling with rates significantly higher than expected decoherence. With the addition of laser fields these can realize tunable-range interacting spin Hamiltonians. We also show how local control of potentials allows isolation of small numbers of ions in a fixed array and can be used to implement high fidelity gates. The use of static trapping fields means that our approach is not limited by power requirements as system size increases, removing a major challenge for scaling which is present in standard radio-frequency traps. Thus the architecture and methods provided here appear to open a path for trapped-ion quantum computing to reach fault-tolerant scale devices.

Figures

Equations

$${\hat{H}}_s$$

$$+\hslash {\omega }_z\left({\hat{a}}_z^{†}{\hat{a}}_z+1/2\right),$$

$${\hat{H}}_{c,ij}$$

$${\Omega }_{\nu ,\mathrm{ex}}^{ij}=\frac{e^2}{4\pi {ϵ}_{0}m{\omega }_{\nu }^{\prime }{R}_{ij,0}^3},$$

$$L=\sum _{j=1}^N\left[\frac{1}{2}m|{\dot{\mathbf{R}}}_j{|}^2+e{\mathbf{A}}_j\cdot {\dot{\mathbf{R}}}_j-e{\Phi }_j\right],$$

$$M\ddot{q}-W\dot{q}+\Phi q=0,$$

$$[{\omega }^{2}M+\omega (-iW)-\Phi ]q_0=0,$$

$$q(t)=\text{Re}\left[\sum _{\lambda =1}^{3N}{r}_{\lambda }{q}_{\lambda }{\mathrm{e}}^{-i\left({\omega }_{\lambda }t+{\delta }_{\lambda }\right)}\right],$$

$$E_{\lambda }$$

$${\hat{a}}_{\lambda }^{†}=\sum _{j=1}^{3N}({\alpha }_{\lambda j}{\hat{p}}_j+{\beta }_{\lambda j}{\hat{q}}_j),$$

$$[{\omega }^{2}M+\omega (-iW)-\Phi ]{\alpha }_{\lambda }=0,$$

$$[{\omega }^2\cdot {\mathbb{I}}_{3N}+\omega (-i{W}^{\prime })-{\Phi }^{\prime }]q_0=0,$$

$$Av=\omega v,$$

$$A=\left[\begin{array}{cc}{\mathbb{O}}_{3N}& {\mathbb{I}}_{3N}\\ {\Phi }^{\prime }& i{W}^{\prime }\end{array}\right].$$

$$\sum _{\lambda =1}^{6N}{\omega }_{\lambda }^2=\text{tr}(A^2)=\text{tr}(2{\Phi }^{\prime }-W^{\prime 2})=-\text{tr}(W^{\prime 2}).$$

$$\sum _{\lambda =1}^{3N}{\omega }_{\lambda }^2$$

$$=N{\omega }_c^2,$$

$${\omega }_{+}^2+{\omega }_{-}^2+{\omega }_z^2={\omega }_c^2.$$

$$\prod _{\lambda =1}^{3N}\left(m{\omega }_{\lambda }^2\right)=|\Phi |.$$

$${\Psi }_{\mathrm{RF}}^{(2)}=\frac{e}{2m{\Omega }_{\text{RF}}^2}{\left[{\Pi }^{(2)}\right]}^2$$

$$\Pi =\frac{{\varphi }_0}{h^2}\left(z^2-\frac{x^2+y^2}{2}\right).$$

$$\Vert {\Psi }_{\text{RF}}^{(2)}\Vert =\frac{\sqrt{3}}{8}|q_z|\cdot \Vert {\Pi }^{(2)}\Vert .$$

$${\hat{H}}_{\text{ODF}}\approx \sum _{j=1}^N{E}_O{\mathbf{k}}_R\cdot {\mathbf{r}}_j\sin ({\varphi }_j-{\mu }_{R}t){\hat{\sigma }}_j^z,$$

$${\hat{H}}_{\text{SPIN}}=\sum _{j{j}^{\prime }}{J}_{j{j}^{\prime }}(t){\hat{\sigma }}_j^z{\hat{\sigma }}_{j^{\prime }}^z,$$

$$J_{j{j}^{\prime }}^0$$

$${\hat{H}}_s$$

$$+\frac{{\hat{p}}_z^2}{2m}+\frac{1}{2}m{\omega }_z^2{\hat{z}}^2$$

$$\hat{x}=\sqrt{\frac{\hslash }{m{\omega }_1}}\left({\hat{a}}_x^{†}+{\hat{a}}_x\right)$$

$$\hat{y}=\sqrt{\frac{\hslash }{m{\omega }_1}}\left({\hat{a}}_y^{†}+{\hat{a}}_y\right)$$

$$\hat{z}=\sqrt{\frac{\hslash }{2m{\omega }_z}}\left({\hat{a}}_z^{†}+{\hat{a}}_z\right)$$

$${\hat{H}}_s$$

$${\kappa }_j=\sum _{k\ne j}\frac{e}{4\pi {ϵ}_0|{\mathbf{R}}_j-{\mathbf{R}}_k|}=k_{e}e\sum _{k\ne j}\frac{1}{|{\mathbf{R}}_{jk}|},$$

$$L_{\text{tot}}=\sum^{N}_{j=1}\bigg{\{}\frac{1}{2}m_{j}|\dot{\mathbf{R}}_{j}|^{2}+e\mathbf{A}_{j}\cdot\dot{\mathbf{R}}_{j}-e\Phi_{j}\bigg{\}},$$

Full text

The study of many-body physics in quantum mechanics is hindered by the inability of classical computing devices to store and manipulate the information required to specify these systems beyond about 50 spins Feynman (1982); Häner and Steiger (2017). Quantum devices, which directly work in Hilbert space, would overcome these limitations and furthermore open up access to other calculations beyond the reach of classical supercomputers Feynman (1982); Georgescu (2014); Schaetz et al. (2013); Blatt and Roos (2012). Trapped atomic ions are among the most successful platforms for exploring these advances. They interact via the long-range Coulomb force, which when combined with laser and microwave fields has allowed high quality quantum logic gates Ballance et al. (2016); Gaebler et al. (2016) as well as Hamiltonian engineering for quantum simulations Zhang et al. (2017). The most successful approach to controlling multiple trapped-ion qubits has involved the use of semi-rigid ion crystals, formed through the competing energy requirements of the global trapping potential and the Coulomb repulsion. Primarily these use 1-dimensional ion chains in radio-frequency (r.f.) ion traps Monz et al. (2011); Zhang et al. (2017). The restriction to one dimension is imposed by the desire to trap at the null of the r.f. field. Penning traps have been used to perform quantum simulations with 2-dimensional (2-d) ion crystals, with the complication of continuous rotation in the crystal plane Britton et al. (2012). In both cases, the intrinsic link between the lattice structure and the oscillation frequencies of the normal modes of oscillation places constraints which limit the range of physics which can be investigated. Furthermore, neither of these approaches is well-suited to scaling these systems up to levels close to a million qubits which are expected to be required for solving relevant problems in quantum chemistry Wecker et al. (2014).

An alternative platform which increases flexibility is the use of micro-fabricated ion traps with electrode structures on length scales close to the inter-ion separation. This would allow access to arbitrary 2-d lattice geometries, as well as providing the possibility to locally tune potentials in order to decouple subsets of ions to facilitate local multi-qubit gates. For r.f. traps, early experiments in this direction have been carried out Wilson et al. (2014a); Mielenz et al. (2016), but a number of significant challenges to scaling arise from the use of radio-frequency potentials. One is that r.f. power dissipation in the electrodes increases with the number of sites (similar to the challenge in scaling optical power which limits optical dipole traps for neutral atoms). It also appears difficult to achieve small site spacings for strong Coulomb coupling in extended arrays of ions Schmied et al. (2009, 2011); Krauth et al. (2015). A further practical challenge is that the operation of such traps relies on the alignment of the microscopic static and r.f. quadrupole potentials, which is challenging to achieve in the presence of stray charges on the electrode surfaces, especially when large numbers of sites are involved Berkeland et al. (1998).

Surface electrode Penning traps offer an alternative platform for trapping ions which utilizes only static fields. Surface-electrode Penning traps have been experimentally realized with both electrons Goldman and Gabrielse (2010) and atomic ions Crick et al. (2010), with the light mass of the former offering the potential for strong site-site couplings Ciaramicoli et al. (2003); Stahl et al. (2005); Goldman and Gabrielse (2010). For electrons, arrays of such traps have been proposed as a coupled system suitable for quantum information Stahl et al. (2005). However electrons lack many of the control techniques available to atomic ions, which have the considerable advantage of being able to utilize advanced laser techniques for cooling, initialization, detection and control Bushev et al. (2008); Hellwig et al. (2010).

In this Letter, we consider the use of arrays of surface Penning traps to realize 2-d lattices of trapped atomic ions for quantum simulation and computation, a setting which appears scalable due to the use of only static fields for trapping. We present a detailed study of the collective motional modes of oscillations for ions trapped in this architecture, which is the crucial element for any quantum information application. We find that the motion of $N$ coupled ions can be mapped onto the same quadratic eigenvalue problem for both the classical and quantum treatments, and use this framework to generalize the single-ion invariance theorem for stable modes of motion Brown and Gabrielse (1982, 1986) to an arbitrary number of ions. We then focus on two cases relevant to quantum simulation and computation in extensible ion lattices. For quantum simulation, we demonstrate that closely spaced arrays can be produced, by obtaining optimal electrode geometries for a range of different lattices. Here Penning traps allow significantly higher ion densities and thus higher coupling than radio-frequency traps for similar experimental constraints. We show that variable range spin-spin interactions can be produced, as well as verifying the ability to access useful ion temperatures through laser cooling. We follow this by studying the feasibility of quantum computing using multi-qubit gates on two nearest-neighbor ions which are embedded in an array, by using local mode decoupling of the ions. We show that a special regime of large zero-point motion available for the magnetron and cyclotron modes allows fast high fidelity logic gates, which could be applied to lattice-based error-correction codes such as the surface and topological color codes Fowler et al. (2012); Bombin and Martin-Delgado (2006).

I Penning microtrap array

The configuration that we propose is shown in figure 1. Electric potentials are applied to trap electrodes laid out on a planar structure such as to form an array of static electric 3-d quadrupoles at a distance $h$ above the trap surface. These, combined with the homogeneous magnetic field, provide three-dimensional confinement of a single ion at each site. Together with the trapping potential the Coulomb interaction defines the equilibrium positions of the ions. For an infinite lattice, the equilibrium positions will align with the centers of the quadrupole potentials created through the periodic arrangement of electrodes. For the finite case, this is not true and should be taken into account while designing the electrodes to generate the required lattice of equilibrium positions. The vibrations of the ions around their equilibrium positions is also coupled by the Coulomb interaction. As we will show in the next section, at the single quantum level this depends as the inverse cube on the distance between ions and the zero-point motion amplitude, and has the strongest effect for two resonant sites. From the trapping perspective, the primary objective is to achieve closely spaced ions. The desirable electrode geometries for infinite lattices will be considered in more detail in section VII, and then applied to quantum simulation problems. A secondary aspect is that the electrodes are microscopic, which allows local control. In this way, ions may be tuned in and out of resonance with the rest. This will be studied in the context of quantum computing in section VIII. In section IX we present a realistic scenario in which these ideas could be implemented in an experiment.

II Dipolar ion-ion coupling

The basis for multi-qubit quantum control and quantum simulations in such a setting is the long-range Coulomb interaction, which couples the oscillations of individual ions at different sites. While for an ion oscillating along a single spatial axis it is simple to think of the oscillating charge distribution as an oscillating dipole potential, and thus that the ions are coupled through something similar to a dipole-dipole coupling, the types of motion exhibited by ions in a Penning trap are more complex. To get a feeling for the nature of this coupling, we first consider for the moment a simplified setting, in which each ion is trapped in a cylindrically symmetric static quadrupole potential $m\omega_{z}^{2}/2(z^{2}-(x^{2}+y^{2})/2)$ for ions of mass $m$ and charge $e$ embedded in a magnetic field of strength $B_{0}$ aligned along the $z$ axis.

At a single site, the Hamiltonian for an ion in such a potential can be written as

[TABLE]

where $\omega_{\pm}=(\omega_{c}\pm\omega_{1})/2$ are the frequencies of the modified cyclotron motion and magnetron motion respectively, with $\omega_{1}=\sqrt{\omega_{c}^{2}-2\omega_{z}^{2}}$ and $\omega_{c}=eB_{0}/m$ the bare cyclotron frequency Brown and Gabrielse (1986). Here $\hat{a}_{\pm}$ and $\hat{a}_{z}$ are the annihilation operators for the corresponding modes. This Hamiltonian differs from a standard 3-dimensional oscillator due to the negative sign in front of the magnetron term. It is also worth noting that both the magnetron and cyclotron motions are two-dimensional, thus the relevant creation and annihilation operators are made up of sums of the form $\hat{a}_{\pm}=\frac{1}{\sqrt{2}}\left(\hat{a}_{x}\pm i\hat{a}_{y}\right)$ where $\hat{a}_{x},\hat{a}_{y}$ refer to motion along $x$ and $y$ respectively (full details of the transformations and definitions used can be found in Appendix A).

Let us now consider two such sites labelled by indexes $i$ and $j$ containing ions with equilibrium positions separated by the vector ${\bf R}_{ij,0}=R_{ij,0}(\sin(\theta_{ij})\cos(\phi_{ij}),\sin(\theta_{ij})\sin(\phi_{ij}),\cos(\theta_{ij}))$ which has magnitude $R_{ij,0}$ and makes an angle $\theta_{ij}$ with the magnetic field. For the current argument, let us work in the approximation that the motion of the ions can be assumed to be a small perturbation which is well-described using a second order expansion of the Coulomb interaction about the equilibrium positions. Moving to a rotating frame with respect to $\hat{H}_{s}$ for the operators at each of the sites and further assuming that the difference frequency between the bare modes is much larger than the respective exchange frequencies for the different modes, we find the Coulomb interaction Hamiltonian

[TABLE]

where $K_{z}=-K_{\pm}=1-3\cos^{2}(\theta_{ij})$. The first term corresponds to hopping of excitations between the sites, while the second gives the modification of the on-site energy due to the static potential of the other ion. The respective coupling strengths for hopping of vibrational quanta from one ion to another are given by the exchange frequencies

[TABLE]

where $\omega^{\prime}_{z}=\omega_{z}$ and $\omega^{\prime}_{+},\omega^{\prime}_{-}=\omega_{1}$.

The coupling Hamiltonian has a dipolar form for all modes. The sign of the coupling for the magnetron and cyclotron modes is inverted with respect to that of the axial motion. For each type of motion the orientation of the effective dipole is along the magnetic field axis. While this is expected for the axial oscillation, for the other modes this is less intuitive, since it corresponds to a direction perpendicular to the plane of oscillation of both the cyclotron and magnetron motions. When the coupling above is generalized to a two-dimensional lattice of sites, the anisotropy of the interactions in the lattice plane depends on the angle $\Theta$ between the magnetic field and the lattice normal. For $\Theta=0^{\circ}$ the interactions are isotropic because $\theta_{ij}=\pi/2$ for all directions within the plane. For $\Theta=90^{\circ}$, $\theta_{ij}$ values vary between [math] and $2\pi$, thus the interactions are anisotropic.

An additional feature of the couplings, which we make use of in section VIII is that the couplings for the modified cyclotron and magnetron motions are proportional to the reciprocal of $\omega_{1}$, rather than the mode frequency itself, a result which stems from the dependence of the zero-point motion of each of these modes. Thus by tuning $\omega_{1}$ to be small, the coupling between two sites can be enhanced. This can be performed by raising $\omega_{z}$ towards $\omega_{c}/\sqrt{2}$ (although the limit $\omega_{z}=\omega_{c}/\sqrt{2}$ is unstable). Unlike a standard mechanical oscillator, this enhancement of the zero-point motion does not require lowering of the mechanical oscillator frequency.

III Operating conditions

Within the restrictions of using symmetric quadrupole potentials, there is nevertheless a range of possible choices of $\omega_{c},\omega_{z}$ which could be used, related to the type of couplings which are desired and the achievable experimental constraints. However these choices are constrained by the fixed dependence of the mode frequencies on these two parameters. As an example which we find interesting from the perspective of engineered tunable-range spin-spin couplings (see section VII.2), we consider the challenge of operating a trap with a fixed achievable $\omega_{z}$ while desiring that the modes have a large enough splitting that three separated bands of modes occur corresponding to coupled magnetron, modified cyclotron and axial frequencies. We find that for an infinite lattice the width of any one band is in the range between 3 and 6 times the exchange frequency for the respective mode. Figure 2 shows the trap frequencies which can be achieved normalized to the axial frequency $\omega_{z}$. To achieve large mode separations, it appears attractive to work in the regime for which $\omega_{+}>\omega_{z}$, which corresponds to $\omega_{c}>3\omega_{z}/2$. However note that this implies that $\omega_{-}<\omega_{z}/2$. The splitting between the $\omega_{\pm}$ modes and the axial modes becomes equal in magnitude when $\omega_{c}=2\omega_{z}$ (see dashed line), for which $\omega_{-}=\omega_{z}(\sqrt{2}-1)/\sqrt{2}\approx 0.29\omega_{z}$. In practice there are good reasons to work with high magnetron frequencies, for example because this would be expected to result in lower heating rates of those modes due to sampling noise at higher frequencies - this may limit the choice of trap parameters.

A second regime of operation could also be used, in which the difference in mode frequencies is significantly below the trap frequencies themselves. For this regime, the axial frequency can be chosen to be the highest frequency of the three, and the modified cyclotron and magnetron modes become closer in frequency. This can be achieved in a narrow regime for which $3\omega_{z}/2>\omega_{c}>\sqrt{2}\omega_{z}$, as indicated by the shaded area in figure 2. As $\omega_{c}$ is reduced within this range, the difference frequency $\omega_{1}$ between the modified cyclotron and magnetron frequencies is reduced. This increases the exchange frequencies for these modes (but also proportionally the heating rate). The maximal separation of modes within this regime is $0.191\omega_{z}$.

IV Normal Mode Analysis

A full analysis of the normal modes of vibration of ions is a prerequisite to understand the methods utilized for laser cooling and implementation of analog quantum simulation. The presence of a magnetic field in Penning traps makes this analysis non-trivial. While the case of a naturally formed two-dimensional ion crystal in a macroscopic Penning trap has been treated before, this involves a transformation to a frame co-rotating with the rigid ion crystal Wang et al. (2013). Here we avoid any such transformation, and derive the normal modes for the general case of ions in an array of micro-Penning traps, assuming only that the equilibrium positions are well defined so that the harmonic approximation for the electric potential can be used and that the magnetic field is uniform over the region explored by each ion. The analysis leads to a generalization of a well known invariance theorem used for single ions in Penning traps Brown and Gabrielse (1982, 1986).

IV.1 Classical Treatment

In the classical regime, the normal mode analysis of a finite system of $N$ trapped ions can be carried out with the help of Lagrangian mechanics. For simplicity, it is assumed that all ions have an identical mass $m$, but the analysis can be generalized to systems containing ions with different masses (see Appendix B). In the frame of reference of the laboratory, with no oscillatory fields present, the Lagrangian of the system is given by

[TABLE]

where $\mathbf{R}_{j}$ denotes the lab coordinates of ion $j$, $\mathbf{A}_{j}=\frac{1}{2}(\mathbf{B}\times\mathbf{R}_{j})$ is the vector potential in the symmetric gauge due to the uniform magnetic field $\mathbf{B}=B_{0}\hat{\bf e}_{z}$, and $\Phi_{j}$ is the electric potential containing a sum of contributions from the trapping potential for ion $j$ and the Coulomb interaction potential experienced by this ion due to all others.

The first step in obtaining normal modes is to solve for the equilibrium positions, which we carry out numerically. The second order term in the series expansion of the system Lagrangian about the equilibrium positions dictates the normal mode dynamics of the system near the stable spatial configuration. The equations of motion for the $3N$-dimensional vector $q=\left[x_{1}\,...\,x_{N}\quad y_{1}\,...\,y_{N}\quad z_{1}\,...\,z_{N}\right]^{T}$ consisting of all the generalized position coordinates can then be deduced as

[TABLE]

where $M$, $W$ and $\Phi$ are $3N\times 3N$ matrices defined as $M=m\cdot\mathbb{I}_{3N}$, $W^{xy}_{jk}=-W^{yx}_{jk}=eB_{0}\delta_{jk}$ and $\Phi^{\mu\nu}_{jk}=\partial q^{\mu}_{j}\partial q^{\nu}_{k}L$. $\Phi$ contains only terms from the static electric potential and the Coulomb interactions. Here, the indices $j,k$ run over the ion numbers $1$ to $N$ while the indices $\mu,\nu$ run over the Cartesian components $x,y,z$. The ‘mass matrix’ $M$ is a real diagonal matrix, $W$ is a real antisymmetric matrix representing the velocity-dependent forces (often referred to as the ‘damping matrix’), while the ‘stiffness matrix’ $\Phi$ is a real symmetric traceless matrix.

To find the normal modes of motion, we substitute the oscillating trial solution $q=q_{0}\mathrm{e}^{-i\omega t}$ which yields the Quadratic Eigenvalue Problem (QEP)

[TABLE]

that can be solved for complex eigenvectors $q_{0}$ and eigenvalues $\omega$. When all eigenvalues are real the motion of all ions is bounded and stable confinement can be achieved. Each of the $3N$ collective normal modes of motion is thus characterized by the eigenpair $\{\omega_{\lambda},q_{\lambda}\}$ and the general solution for the motion can be expressed as

[TABLE]

with the amplitude $r_{\lambda}$ and phase $\delta_{\lambda}$ for each mode $\lambda$ determined by the initial conditions.

It is important to note that the total energy contained in each mode

[TABLE]

is not trivially positive, unlike the case of Paul traps. Typically we will observe $N$ modes dominated by motion along the axial direction and it is convenient to continue calling these modes axial modes in the context of the $N$ ion array of Penning traps. Each of the axial modes has a positive total mode energy. Similarly there are $2N$ radial modes out of which $N$ have each a positive mode energy and $N$ have each a negative mode energy. We will denote the radial modes with positive sign as cyclotron modes and the ones with negative sign as magnetron modes.

IV.2 Quantum Mechanical Treatment

The solution for the normal modes in the quantum regime involves the formulation of the Hamiltonian in terms of the canonical position and momentum operators $\hat{q}_{j}$ and $\hat{p}_{j}$, and then forming the phonon creation and annihilation operators, $\hat{a}^{\dagger}_{\lambda}$ and $\hat{a}_{\lambda}$, for each mode $\lambda$ as linear combinations of these operators,

[TABLE]

where $\alpha$ and $\beta$ are complex coefficients. The objective is to find these coefficients which allow us to diagonalize the Hamiltonian for a stable system in the second quantized form $\hat{H}=\sum^{3N}_{\lambda=1}\hbar\omega_{\lambda}(\hat{a}^{\dagger}_{\lambda}\hat{a}_{\lambda}+\frac{1}{2})$ with the phonon operators following the standard commutation relations, $[\hat{a}^{\dagger}_{\lambda},\hat{a}^{\dagger}_{\lambda^{\prime}}]=0$, $[\hat{a}_{\lambda},\hat{a}_{\lambda^{\prime}}]=0$, $[\hat{a}_{\lambda},\hat{a}^{\dagger}_{\lambda^{\prime}}]=\delta_{\lambda\lambda^{\prime}}$. Combining these requirements, we find that the $3N$-dimensional vectors $\alpha_{\lambda}$ for each mode $\lambda$ satisfy the same QEP we had to solve in the classical analysis

[TABLE]

and we have the relation $\beta_{\lambda}=i\omega_{\lambda}M\alpha_{\lambda}+\frac{1}{2}W\alpha_{\lambda}$. These vectors can then be normalized such that the condition $[\hat{a}_{\lambda},\hat{a}^{\dagger}_{\lambda}]=1$ is fulfilled.

We note that the treatment discussed in this section encompasses both r.f. and Penning traps, and the normal modes for the former under the pseudopotential approximation can be retrieved by choosing the magnetic field strength as zero and adding a suitable pseudopotential term to $\Phi$. In this case $\Phi$ is not trace zero.

V Generalized Invariance Theorem

A real trap is imperfect and can suffer from misalignments between the magnetic field and the confining axis of the quadrupole potential, or the trap potential differing from the desired precise form. Including these imperfections in the matrices $W$ and $\Phi$ in the QEP encountered in the classical analysis and dividing the equation by $m$, we get

[TABLE]

where we define the matrices $W^{\prime}=W/m$ and $\Phi^{\prime}=\Phi/m$. This new QEP can then be linearized by mapping it on to a standard eigenvalue problem while increasing the dimensionality by a factor of two so that we arrive at

[TABLE]

with $6N$-dimensional eigenvectors $v=\begin{bmatrix}q_{0}&\omega q_{0}\end{bmatrix}^{T}$ and $6N$ eigenvalues $\omega$ belonging to the $6N\times 6N$ matrix $A$

[TABLE]

Since $A^{2}v=\omega^{2}v$ and the sum of eigenvalues of a matrix is equal to its trace,

[TABLE]

This result will hold for any potential terms added to $\Phi^{\prime}$ which exist in free space, since these must be traceless in order to satisfy Laplace’s equation (it would not hold for a pseudopotential). While for a single ion only the trapping potentials are present, in the case of a multi-ion system the Coulomb interactions are also contained in this term. Noting that the frequencies come in pairs of positive-negative values in the stable regime we can express this sum in terms of the $3N$ positive frequencies,

[TABLE]

where $\omega_{c}=eB_{0}/m$ is the bare cyclotron frequency. This relation between the strength of the magnetic field quantified and the normal mode frequencies of a stable $N$-ion system represents a non-trivial generalization of the well known Brown-Gabrielse invariance theorem for a single ion in a Penning trap Brown and Gabrielse (1982, 1986),

[TABLE]

A complementary result relating the product of the normal mode frequencies to the curvature of the total electric potential can be derived by taking the determinant of the matrix $A$, giving

[TABLE]

The results (15) and (17) further generalize to systems containing ions of different masses (see Appendix D for more details). The single-ion invariance theorem is widely used in precision measurement Gabrielse (2009), and our generalization can be applied to precision mass measurements of multi-ion crystals Gutiérrez et al. (2019).

VI Laser Cooling

Doppler cooling is more complicated in Penning traps as compared to radio-frequency traps due to the fact that the magnetron modes have a negative total energy. As a consequence the cooling requirements of the magnetron modes are incompatible with those of the axial and cyclotron modes, and no combination of uniform beams can cool all modes of motion simultaneously Itano and Wineland (1982). One way to combat this is to use a non-uniform beam with intensity gradient, but the final temperatures reached for both kinds of radial modes are greater than one would expect from the standard Doppler cooling limit, and the range of motional frequencies that allow for cooling all three kinds of modes is also restricted Itano and Wineland (1982).

An alternative solution which has been realized experimentally Powell et al. (2002) is to couple the cyclotron and magnetron modes by applying a weak quadrupolar electric field $\phi_{\rm ax}(x^{2}-y^{2})$ oscillating at the bare cyclotron frequency, in a technique known as axialization Thompson and Papadimitriou (2000). A red-detuned uniform-intensity Doppler cooling laser beam can then simultaneously cool all modes. With the axialization drive the system no longer consists solely of electrostatic fields but since the amplitude of such a drive is much lower than for the r.f. drive required for Paul traps, the deleterious effects of micromotion are accordingly much smaller. Moreover, this technique is required only during the laser cooling process and works efficiently at all trap frequencies, allowing trapping in regimes not accessible through the use of just inhomogeneous beams.

The derivation of the rate equations of the mode amplitudes of ions due to Doppler cooling in Penning traps in the presence of axialization has no simple analytical solution. Instead, we perform numerical simulations on small numbers of ions in which we numerically integrate the equations of motion of the trapped ion including the axialization potential as well as stochastic momentum kicks which occur with a probability which depends on the Doppler shift between the laser and the ion resonance to simulate photon scattering events. By running the simulation a large number of times, the average classical amplitudes of each mode can be found, which can then be converted to mean quantum phonon occupation numbers. Results of such simulations for laser cooling of a six 9Be+ ion honeycomb lattice with lattice constant $15$ µm with the confining axes tilted at $\Theta=20^{\circ}$ with respect to the radial plane are shown in figure 3. Here $B_{0}$ = $2.5$ T and $\omega_{z}=2\pi\cdot 2.1$ MHz. The uniform laser beam is oriented parallel to the electrode plane so that $\mathbf{k_{L}}=\cos\Theta\hat{\bf e}_{x}+\sin\Theta\hat{\bf e}_{z}$. The axialization voltage here is $\phi_{\rm ax}=0.03\phi_{0}$. The initial quantum numbers for each mode are chosen within a range close to $10^{4}$ quanta. The results show that modes can be cooled with time constants in the range of 0.1-0.2 ms, with final mode occupations in the range of 10-30 quanta. This would be expected for Doppler cooling at these trap frequencies. For coupled ion arrays, it is important that the axialization drive mixes all magnetron modes with the modified cyclotron modes. This implies that the modulation strength must be greater than the width of the relevant spread of frequencies of each set of modes.

VII Quantum simulation in fixed lattices

One of the primary uses of multi-ion control in recent years has been for quantum simulation of lattice models. This involves the use of always-on couplings between either spins or motional degrees of freedom of the different ions, realizing a representation of a Hamiltonian of interest elsewhere in physics. Although proposals and experimental demonstrations exist for spin-boson and pure bosonic systems Porras et al. (2008); Barreiro et al. (2011); Bermudez et al. (2011), a particular focus of recent work has involved Ising spin models implemented on the internal degrees of freedom Britton et al. (2012); Zhang et al. (2017). In the following section, we describe how Penning trap arrays could be used for such studies, including optimal electrode layouts, the relevant features of the normal modes, and the implementation of tuneable range Ising spin interactions. These ingredients form a hierarchy; lattice, modes and spin Hamiltonians which are largely shared with other types of quantum simulations which might be of interest. A strength of the systems considered here is that the lattice can be designed through the electrode layout.

VII.1 Optimal Electrode Geometries

From equation 3 above, it is clear that the $1/{R^{3}_{ij,0}}$ nature of effective couplings between the ions favors forming closely-spaced ion arrays. Although this can be achieved by scaling the size of the whole trapping structure, this results in a reduction of the ion-electrode distance which is undesirable due to the expected increase in ion motional heating Deslauriers et al. (2006); Brownnutt et al. (2015), and the increased chance that stray scattered light from the optical control fields used for cooling and engineered spin-spin interactions induces charging of the electrode surfaces resulting in stray electric fields. For operation of the system, it is also desirable to work with trap frequencies which are high enough to avoid common noise sources in the laboratory, which reduces heating and facilitates laser cooling. For a given electrode structure, the motional frequencies can be increased with a corresponding increase in the applied electrode voltages. However at some point this will be limited by voltage breakdown, and therefore it is beneficial to search for optimal electrode layouts for achieving closely spaced ion traps while retaining high trap frequencies of the individual micro-traps. We consider here the experimental feasibility to generate such surface-electrode trap layouts with high motional coupling between ions in micro-Penning trap arrays for a given applied voltage. Our focus lies in particular on single layer surface-electrode traps as they offer an open planar structure which facilitates optical access. We note that approaches with two planes of electrodes facing each other might allow improved conditions with regards to spin-spin couplings, but these seem to be more technically challenging Krauth et al. (2015).

Previous work has described methods for obtaining surface-electrode geometries which maximize the achievable curvature of the pseudo-potential in arrays of r.f. traps for a given trapping site density with distance from the electrode surface $h$ Schmied et al. (2009). The problem reduces to maximizing the quadrupole strength which can be produced at the array of sites, which is the identical problem for the Penning trap arrays. However in the case of r.f. traps this quadrupole potential must be converted into a ponderomotive pseudopotential while maintaining conditions suitable for stable motion, while the Penning trap frequencies are directly dependent on the static quadrupole. The advantage this gives can be evaluated by considering the effect of modulating a static potential $\Pi$ with curvature tensor $\Pi^{(2)}\equiv\partial_{\mu}\partial_{\nu}\Pi$, with $\nu,\mu=x,y,z$ at a radio-frequency $\Omega_{\rm RF}$, creating $\Pi\cos(\Omega_{\rm RF}t)$. In the pseudopotential approximation Wineland et al. (1998), the curvature tensor of the pseudopotential for an ion of mass $m$ is then

[TABLE]

at any trap center. For simplicity we focus on a cylindrically symmetric trap potential defined as

[TABLE]

To compare the strength achievable for the r.f. vs Penning trap, we take the Frobenius norm of the curvature tensor, finding that the magnitudes of the two curvature tensors can be related using the Mathieu parameter $q_{z}=-4e\phi_{0}/\left(m\Omega^{2}_{\text{RF}}h^{2}\right)$ as

[TABLE]

Thus the curvature of the pseudopotential is weaker than that of the corresponding static potential by a factor $\sqrt{3}\lvert q_{z}\rvert/8$. For surface-electrode r.f. traps, stability becomes a concern for $\lvert q_{z}\rvert\gtrsim 0.3$, corresponding to a reduction factor of around 15. In this case, for a given voltage which is applied to the electrodes, the trap frequency in the r.f. trap is reduced relative to a Penning trap with the same geometry by a factor $\sqrt{15}$.

The discussion above makes it clear that the optimization of electrode structures for Penning trap arrays is identical to the case of radio-frequency traps, and thus produces identical electrode geometries Schmied et al. (2009, 2011). Similar to the earlier work, we define a dimensionless curvature $\kappa=\lVert\Phi^{(2)}\rVert h^{2}/V$, where $V$ is a fixed voltage applied to part of the electrode plane (the rest being set to zero), and $h$ is taken to be the distance from the centre of the quadrupole to the nearest trap surface. For the symmetric potential the component of the dimensionless curvature aligned with the magnetic field is $\hat{\kappa}_{z}=\hat{\kappa}2^{2/3}$. We then optimize the electrode geometry. Figure 4 shows the the dimensionless curvature of the trap potential achievable for different infinite lattices as a function of the ratio of trap height to inter-ion distance, with the magnetic field directed perpendicular to the plane of the electrodes. The values given are normalized to the dimensionless curvature $\hat{\kappa}\approx 0.473$ achievable for an optimal surface-electrode point trap Schmied et al. (2009).

Typical electrode structures in quantum information experiments can withstand differences in voltages on neighboring electrodes sufficient to achieve $V=300$ V Sterling et al. (2013). With $h=30$ µm, $\kappa=10^{-4}$ then allows to achieve $\omega_{z}=2\pi\cdot 2.1$ MHz for beryllium ions. Figure 4 indicates that this would allow ions spaced by between $26$ µm and 15 µm, with the former corresponding to the triangular lattice and the latter to the honeycomb and Kagome lattices. The resulting exchange frequencies are between 11 kHz $<\Omega_{\rm ex,z}/(2\pi)<55$ kHz. This is far above heating rates and frequency drift rates observed in traps of a similar size, thus high quality coherent exchange would be expected Brown et al. (2011).

The discussion above considers trapping potentials with the confinement direction (and magnetic field) out of the plane of the electrodes. The introduction of a tilt between the electrode plane normal and the magnetic field modifies the geometries of the optimized electrodes and the values of $\kappa$ which can be achieved. Figure 5 shows the value of $h/d$ for which a dimensionless curvature $\kappa/\hat{\kappa}=10^{-4}$ is achieved as a function of the angle $\Theta$ between the magnetic field and the normal to the electrode plane. We see that $\Theta=90^{\circ}$ produces the highest curvatures for all geometries (and thus requires the smallest $h/d$ to get to $\kappa/\hat{\kappa}=10^{-4}$), with an additional maximum for $\Theta=0^{\circ}$. The former allows laser cooling with laser beams in the plane of the electrode surface, the latter does not - previous work in radio-frequency traps indicates that $\Theta>8^{\circ}$ is necessary for robust laser cooling Lindenfelser et al. (2017). The change in achievable $h/d$ for a fixed $\kappa$ between $\Theta=0^{\circ}$ and $\Theta=8^{\circ}$ is small for the triangular and square lattices, but considerable in the case of the Kagome and honeycomb lattices.

VII.2 Spin-spin interactions for the Ising Model

We now examine the possibility of implementing spin-spin couplings of the form relevant to studying models such as the transverse Ising model Ising (1925), which has been a common target of quantum simulation experiments using trapped ions. An effective spin-spin interaction can be generated based on standard techniques developed in the trapped-ion quantum computing community to implement multi-qubit gates Leibfried et al. (2003); Mølmer and Sørensen (1999); Roos et al. (2008); Kim et al. (2009). These methods rely on the application of forces that depend on the internal (pseudo)spin-state of the ions. For instance, two laser beams off-resonant with respect to the internal transition and with a frequency difference $\mu_{R}$ and wave-vector difference $\mathbf{k}_{R}$ between each other create a traveling-wave interference pattern at the ions. Each ion experiences a state-dependent optical dipole force (ODF) that oscillates at the frequency $\mu_{R}$. To simplify the algebra, we assume that the two relevant states of the ions are spin-half ground states with no hyperfine structure Home et al. (2006). In this case, it is possible to generate an ODF that is equal in magnitude but opposite in sign for the two internal states which can be considered as eigenstates of the Pauli operator, $\hat{\sigma}^{z}$. For small coherent displacements of the ions from their equilibrium positions we can use the Lamb-Dicke approximation and keep only resonant terms, resulting in the ODF interaction Hamiltonian

[TABLE]

where $E_{O}$ depends on the laser beam properties as well as the matrix elements of the internal transition of the ions, and the phase at the ion location is given by $\phi_{j}=\mathbf{k}_{R}\cdot\mathbf{R}_{j0}$. A similar Hamiltonian can be achieved in a rotated spin basis ($\hat{\sigma}^{x},\hat{\sigma}^{y}$) by driving both red and blue motional sidebands of the spin-flip transition simultaneously Roos et al. (2008). Given the periodic arrangement of ions, we can ensure that this phase is the same for all ions using well-chosen laser beam configurations. To simplify the following discussion, we assume that this condition is met (see Appendix E for a more general treatment).

The time evolution operator associated with $\hat{H}_{\text{ODF}}$ can be calculated by carrying out a Magnus expansion in the interaction picture, and for the given ODF interaction this yields two terms. The first term describes periodic spin-motion entanglement generated by the ODF. Quantum simulation experiments typically work in the regime in which this entanglement is negligible and can therefore be adiabatically eliminated Blatt and Roos (2012). The second term describes an effective Ising-like spin Hamiltonian

[TABLE]

with the static part of the spin-spin interactions $J_{jj^{\prime}}(t)$ given by

[TABLE]

Here we have defined the Lamb-Dicke parameters in a slightly unconventional manner as $\eta_{\lambda j}=\sum_{\nu=x,y,z}k^{\nu}_{R}\rho_{\lambda 0}\gamma_{\lambda j\nu}$, with $\rho_{\lambda 0}$ being the spread of the zero-point wavefunction of mode $\lambda$ and $\gamma$ is the corresponding eigenvector normalized to 1. For an r.f. trap this definition reduces to the standard form of the Lamb-Dicke parameter Home et al. (2011).

An interesting aspect of the simulation of the Ising model using trapped ions is the possibility to engineer spin-coupling terms which follow a power-law scaling with the inter-ion separation $\lvert J_{jj^{\prime}}\rvert\propto 1/\lvert R_{jj^{\prime}0}\rvert^{a}$, with $a$ dictated by the experimental detunings Britton et al. (2012). However, such tunable-range interactions are possible only with certain mode structure, in which the center-of-mass (COM) mode has the highest (lowest) frequency in a given band, and the state-dependent force is tuned to a higher (lower) frequency than this mode. In previous experiments with bulk crystals this situation is naturally satisfied, whereas it is not always satisfied in the trap arrays which we consider in this paper. In the following, we trace the importance of the normal mode structure in the determination of the effective spin-spin interactions that can be engineered for a given system of ions. We take as an example a honeycomb lattice of 204 ions, with the nearest-neighbor separation $d=15$ µm. Here $B_{0}$ = $2.5$ T and $\omega_{z}=2\pi\cdot 2.1$ MHz. We use an ODF interaction strength corresponding to a Rabi frequency of $E_{O}/\hbar=2\pi\cdot 300$ kHz, which is a level similar to that used in previous experiments Zhang et al. (2017).

We first consider the case when the magnetic field is aligned normal to the plane in which the ions lie. With such an orientation the axial motion is decoupled from the radial motion and the COM mode sits at the highest frequency in the axial band. With a wave-vector difference $\mathbf{k}_{R}=k_{R}\hat{\bf e}_{z}$, only the axial modes are excited. By tuning the ODF to the blue of the axial branch variable range spin-spin couplings can be engineered with the range of interaction decreasing from infinite range ($a=0$) to dipole-dipole type ($a=3$) as $(\mu_{R}-\omega_{z})$ is increased. Since all coupling terms are positive, this allows to simulate an antiferromagnetic Ising Hamiltonian. Experiments carried so far using both r.f. traps (for eg. Ref. Zhang et al. (2017)) and Penning traps (for eg. Ref. Britton et al. (2012)) are based on this simplification. Conversely, a tunable-range ferromagnetic Ising model can be simulated by aligning $\mathbf{k}_{R}$ along the radial plane and tuning the ODF to the red (blue) of the cyclotron (magnetron) branches. This is possible since the cyclotron (magnetron) COM mode is the lowest (highest) in its branch. The coupling terms for the magnetron branch are negative since the magnetron motion represents an inverted harmonic oscillator.

One of the challenges for realizing a setup with the magnetic field normal to the plane is that it is difficult to cool the axial motion or generate an ODF along the axis using lasers directed parallel to the surface of the chip, or equivalently, the plane where the ions sit. While these problems can be countered by using in-chip waveguides Mehta et al. (2016), we rather attempt to see how well the behavior described above holds when the magnetic field is tilted at an angle $\Theta$ with respect to the normal of the lattice plane. The crucial factor here is the position of the COM mode within the branch being excited by the ODF. For the trap settings and lattice considered here, the COM mode is on the edge of the band for $\Theta\lesssim 43^{\circ}$. Figure 6 shows the normal mode spectrum for the case of $\Theta=20^{\circ}$, which does satisfy this condition. Spin-spin coupling terms generated from this lattice through tuning the ODF outside the frequency spectrum of the cyclotron and axial branches are plotted in figure 7.

For the case of $\Theta=90^{\circ}$, that is, when the magnetic field is along the plane, a suitably oriented in-plane laser beam can cool all motional modes. However all three COM modes lie away from the extrema of their respective branches, making it hard to implement variable range spin-spin interactions by tuning $\mu_{R}$ alone. With such a mode structure, detuning $\mu_{R}$ in either direction from any of the three COM modes does not reveal a well-defined power law decay and the coupling terms have different signs depending on the angle of the inter-ion vector. A histogram plot in figure 8 shows this behavior. The same is expected more generally whenever $\mu_{R}$ lies in the middle of the branch of modes used. More complicated methods involving multi-frequency laser beams could make it possible to emulate a tunable-range Ising model with such geometric arrangement of micro-traps Korenblit et al. (2012), but we do not consider this here. At large detunings of the ODF from a given branch, dipole-dipole couplings can be realized with a distance scaling $\lvert J_{jj^{\prime}}\rvert\propto 1/\lvert R_{jj^{\prime}0}\rvert^{3}$ but the sign of the coupling term between any given pair of spins will be determined mostly by their relative phase in the mode closest to the chosen value of $\mu_{R}$. This frustration in sign might allow for the study of disordered spin dynamics in poorly understood systems such as quantum spin glasses Edwards and Anderson (1975). The same behavior could also be effected by tuning within any phonon branch although the emergence of any power law scaling is not expected, except for the case when the ODF is tuned close to the COM mode, leading to infinite range behavior ($a=0$).

While the statements above give a general discussion, they do not include the tuning of the angle between the projection of the B-field into the plane and the lattice symmetry axes. Here it is probable that special cases arise with interesting features. This is an area of future study.

VIII Quantum computation with local modes

Many of the most promising approaches to quantum error-correction also make use of extended two-dimensional lattices of qubits. These include both the surface code Fowler et al. (2012) and the topological color codes Bombin and Martin-Delgado (2006). For fault-tolerance errors must be local and must have low rates, thus it is desirable that gates between ions involve only the chosen ions and do not require precise control of the complete array. For this reason it is desirable to decouple these chosen ions from the rest of the array. In our architecture this can be achieved by local tuning of the trap potential of the ions in question.

As a test case, we consider a lattice of 90 beryllium ions on a square lattice spaced by 30 µm with a quadrupole potential corresponding to $\omega_{z}=2\pi\cdot 2.55$ MHz (see figure 9 i). We choose a magnetic field of 2.2 T which lies in plane, for which $\omega_{c}/(2\pi)=3.75$ MHz. For an ion in a single isolated potential, the modified-cyclotron and magnetron modes are then at $2.39$ MHz and $1.36$ MHz respectively. In order to perform a gate locally between the central two ions, we tune the curvature of the potential for these two ions such that the axial oscillation frequency for a single isolated ion is $(\omega_{z}+\Delta\omega_{z})=2\pi\cdot 2.63$ MHz, where $\Delta\omega_{z}=2\pi\cdot 80$ kHz. The resulting mode spectrum is shown in figure 9. The choice of curvature is special, because the potential it produces meets conditions for which the resulting modes have a relatively small separation between the uncoupled magnetron and modified-cyclotron modes (a contribution to which is made by the potential of the neighboring ions). Since the zero-point motion of these modes scales as $1/\Delta_{\alpha}$ with $\Delta_{\alpha}=\omega_{\alpha}-\omega_{c}/2$, such a choice enhances the zero-point motion of the modified cyclotron and magnetron modes. In addition, these modes are relatively far-detuned from the modes of the rest of the lattice. For the chosen curvature, the coupled modes closest to $\omega_{c}/2$ are those in which the selected ions move out of phase (often called stretch modes), which have frequencies for which $\Delta_{\rm s,+}=2\pi\cdot 60.2$ kHz and $\Delta_{\rm s,-}=2\pi\cdot-60.50$ kHz and zero-point motion for the chosen ions of $96.4$ nm. The next closest are the center-of-mass (COM) modes of the selected ions, with $\Delta_{\rm c,+}=2\pi\cdot 179.3$ kHz and $\Delta_{\rm c,-}=2\pi\cdot-179.6$ kHz with corresponding zero-point motion of $55.9$ nm. These then form an isolated set of modes on which a multi-qubit gate can be performed. Similar to the approach taken in section VII, we consider a geometric phase gate Leibfried et al. (2003); Home et al. (2006) which uses an oscillating state-dependent force such as can be produced with a traveling standing wave or a magnetic field gradient, and make the simplification that the force at any given point is equal in magnitude and opposite in sign for the two eigenstates of $\hat{\sigma}_{z}$. The Hamiltonian is then the one found in equation 21, and we assume that the experiment can be arranged such that the phase of the force is the same at each ion. To perform a gate, the frequency $\mu_{R}$ must be chosen. For the mode detunings above, an attractive possibility is to use an oscillating force at $\mu_{R}=\omega_{c}/2$, which drives the two stretch modes almost equally, but with the opposite detuning. While for a Paul trap this would result in the phases due to each mode cancelling out, for the Penning trap the contributions of both modes add, because the lower frequency mode of the pair is a magnetron mode and thus has a negative frequency. If only these two modes were included in the gate, it could be performed in a time of $t_{g}=2\pi/\Delta_{s,+}$ by using a Rabi frequency $\eta_{s}\Omega=\Delta_{\rm s,+}/\sqrt{2}$. In practice the contributions of the COM modes subtract from this effect, and the additional bulk modes also contribute. For a Lamb-Dicke parameter of $\eta\sim 0.17$ for the stretch modes, simulations involving all modes show that a Rabi frequency of $E_{O}/\hbar=2\pi\cdot 300$ kHz could be used to perform a gate which would produce a Bell-state with fidelity of $F>0.9998$ in 16 µs. The large zero-point motion means that a Raman beam pair with a small difference wavevector is required to operate within the Lamb-Dicke regime, which is desirable for insensitivity of gate fidelity to initial ion temperature Leibfried et al. (2003). For beryllium this would require an angle of $\theta_{R}=\pi/36$ for the beryllium wavelength of $313$ nm. Alternatively a magnetic field gradient of $\sim 19$ T/m (lower than that realized in recent experiments Khromova et al. (2012)) would provide a gate with the same speed. Note that the motional mode parameters used in this analysis were chosen to satisfy a close-to-integer ratio between $\Delta_{\rm c}$ and $\Delta_{\rm s}$, such that both the local stretch and COM modes are dis-entangled from the internal states at the end of the gate Leibfried et al. (2003).

These results indicate that local changes to the potential, combined with individual optical or microwave addressing of the ions could be used to realize quantum computing in the proposed architecture. The enhanced zero-point motion used here is particularly appealing in the context of magnetic field gradient gates, which struggle to achieve high gate speeds in Paul trap settings due to the challenge of producing high field gradients Ospelkaus et al. (2011); Harty et al. (2016a); Weidt et al. (2016); Khromova et al. (2012). In the presence of high Rabi frequencies, faster gate speeds should be possible using multi-pulse techniques following methods demonstrated in Paul traps Steane et al. (2014); Schäfer et al. (2018). For error-correction the need for regular detection of ancilla ions poses challenges with regards to measurement cross-talk, which might require the use of selective electron shelving Wineland et al. (1998). For parallelizing error-correction codes, it would be necessary to select multiple pairs of ions at different points in the lattice, and perform gates on each of these simultaneously. Here the $1/d^{3}$ nature of the Coulomb mode coupling is advantageous.

IX Implementation example

As an example of a possible implementation of a two-dimensional Penning trap array suitable for quantum simulation we consider a honeycomb lattice of 62 sites. The surface-electrode pattern required for such a layout with nearest-neighbor spacing of 30 µm is shown in figure 10. Applying a voltage of 135 V to the DC electrodes results in an axial trap frequency of $\omega_{z}=2\pi\cdot 2.1~{}\text{MHz}$ for beryllium ions. In a global magnetic field angled at $\Theta=20^{\circ}$ from the plane normal with a magnitude of $B_{0}=2.5$ T, we get a reduced cyclotron frequency of $\omega_{+}=2\pi\cdot 3.73~{}\text{MHz}$ and a magnetron frequency of $\omega_{-}=2\pi\cdot 0.54~{}\text{MHz}$. In this setup, effective laser cooling requires an axialization field, which is generated by applying a drive to the r.f. electrodes at a frequency of $\omega_{+}+\omega_{-}=2\pi\cdot 4.27~{}\text{MHz}$ and an amplitude of 4 V.

Working inside the bore of a several tesla superconducting magnetic presents challenges for the delivery of the optical power required for cooling, detection, and single- and multi-qubit control operations. One approach which has been used in multiple experiments is to direct beams down the magnet bore, and then re-direct them using mirrors through a vacuum system housed in the bore Biercuk et al. (2009a); Mavadia et al. (2013); Ball et al. (2019). For the purposes of working with ions in an array of surface traps with the ions at a height of 45 µm from the surface, we envision beams with a diameter of 2 mm directed down the bore of the magnet, which are then reflected towards the trap with a mirror, and then focused to a beam waist of 10 µm using a 100 mm focal length cylindrical lens. This configuration results in an intensity variation of only 2.5% over the trap area, where 100 µW of 313 nm cooling light would have an intensity 8 times that of the saturation intensity, which would be sufficient for laser cooling. Raman beams with a power of 4.3 mW can be directed parallel to the trap surface with the difference wave-vector chosen to address either the axial or the radial modes of motion. This results in a one-photon Rabi frequency of $\Omega=2\pi\cdot 300~{}\text{MHz}$, and for a detuning of 150 GHz a two-photon Rabi frequency of $\Omega_{2}=2\pi\cdot 300~{}\text{kHz}$.

The final element required for running such a setup is the imaging system. Here we envision using a Schwarzschild objective that is diffraction limited over a 130 µm2 region, with an effective focal length of 9.5 mm, a working distance of 18 mm, and a collection efficiency of 4%, similar to what we have developed for previous experiments Leupold (2015). Light collected from this objective is directed back out of the magnet bore with a mirror where it can be imaged onto a suitable detector. Taking into account other losses and the quantum efficiency of the detector, we would expect a detection efficiency of 1%, or 50 counts in 100 µs.

It is worth considering the power dissipation due to the axialization field. For a surface electrode trap fabricated using a Complementary-Metal-on-Silicon (CMOS) process, each ion’s trap site (with electrode area of 1200 µm2) would be expected to have a capacitance of 0.01 pF. Assuming a resistance of $R=4$ ohms of the wiring to the trap chip this means that each trap site would dissipate $P_{\text{RF}}=\frac{1}{2}RC^{2}V_{\text{RF}}^{2}\Omega_{\text{RF}}^{2}=2.3~{}\text{pW}$ due to the axialization. Such a treatment makes clear the power-dissipation advantage of the Penning trap approach. For a Paul trap of similar dimensions and trap frequencies the required r.f. voltage is $V_{\text{RF}}=100\text{\,V}$ at a drive frequency $\Omega_{\text{RF}}=2\pi\times 163~{}\text{MHz}$ Wilson et al. (2014b). This would dissipate a power of $P_{\text{RF}}=2~{}\text{\textmu W}$, a factor of $10^{6}$ more.

X Scaling

Arrays of microfabricated Penning traps provide a new perspective towards scaling, because they do not require applying r.f. drives to the ion trap chip beyond the modest frequencies and voltages needed for axialization. The reduction in power-dissipation becomes even more critical for a larger number of sites, for instance an array similar to that considered above with thousands of sites (as might be required for quantum computation) would require dissipating several hundred watts in a similarly sized r.f. trap, which is most likely not feasible. By contrast, the Penning trap approach dissipates several hundred microwatts due to axialization fields, which only have to be applied during the reduced fraction of time devoted to Doppler cooling.

Penning traps also reduce the number of degrees of freedom which must be controlled at any one site relative to radio-frequency traps. Stray electric fields in Paul traps lead to misalignment of the radio-frequency and static quadrupole potentials, resulting in undesirable micromotion which affects the interaction of the ions with laser fields. In Penning traps, these move the center of the trap, resulting in shifts of the motional frequencies if the potential is anharmonic, but do not produce any other undesirable effect. Recent evidence suggests that heating may be linked to processes driven by the radio-frequency drive in Paul traps Hite et al. (2017), which gives hope that anomalous heating could be reduced in the systems proposed here.

While the prospects for scaling look attractive relative to current approaches, many caveats remain. Although co-wiring of thousands of electrodes is well within the capabilities of fabrication based on commercial CMOS processes, and trap chips fabricated using these methods have been successfully used to trap ions Mehta et al. (2014), much work on fabrication and operation remains in order to realize large-scale arrays. In common with all ion trapping schemes requiring local control for quantum computing, a considerable remaining challenge is optical delivery of focused laser beams, in particular where local fields for individual sites are required (this constraint is not met in the approach considered in section IX above, which uses global laser beams). For truly large-scale quantum computing systems the integration of optics Mehta et al. (2016, 2020); Niffenegger et al. (2020) into the ion trap chips seems to be essential. The need for high-intensity laser beam delivery might be mitigated using oscillating magnetic fields delivered directly from the ion trap chip, or using static magnetic field gradients Srinivas et al. (2019); Harty et al. (2016b).

As opposed to Paul traps, the primary difficulty of realizing high-quality qubit control in Penning traps is the presence of a high magnetic field, fluctuations in which pose a limit to spin-state coherence. Stability on the level of a part-per-billion is available through the use of superconducting magnets, which limits current experiments to coherence times of 50 ms Britton et al. (2016). The effect on the stability of motional modes is negligible at this level. One approach to protecting spin qubits in a quantum computation setting would be to use decoherence-free-subspaces, which protect against the expected homogeneous field fluctuations Häffner et al. (2005). Alternative schemes might utilize dynamical decoupling or employ nuclear transitions (which are less sensitive to field fluctuations) for memory Biercuk et al. (2009b).

XI Conclusions

This study establishes the possibility of using ions in Penning trap arrays for scalable many-body quantum simulations and quantum computation. While we consider here simple settings for motional modes and the possibility to realize tunable-range spin-spin interactions, the flexibility of local control of trapping potentials means that couplings could be used to access a wide range of possibilities which have been previously discussed in the context of other systems, including but not limited to spin-boson systems Porras et al. (2008), dissipative simulations Barreiro et al. (2011), and engineered topology Bermudez et al. (2011). Although quantum computing seems feasible using static arrays with selected ions tuned into local resonance, this is only one way of scaling trapped-ion quantum information. Breaking the lattice down into smaller units would allow smaller ion separations to be achieved Mielenz et al. (2016). Utilizing this would require some level of transport of ions, which could be performed by moving the electric quadrupole positions dynamically Wineland et al. (1998); Kielpinski et al. (2002); Hellwig et al. (2010); Crick et al. (2010). Here the Penning microtrap array holds the considerable advantage that there is no need to separate regions for quantum gates from specialized junction regions where 2-dimensional transport occurs Blakestad et al. (2009); Shu et al. (2014). Since the homogeneous magnetic field supplies 3-dimensional confinement anywhere that a static quadrupole can be placed, re-organization of the potential landscape in 3 dimensions would allow 3-dimensional movement of ions at any point above the trap surface. Thus Penning microtrap arrays appear to remove multiple constraints on scaling trapped-ion quantum computing, paving the way to useful quantum computers.

Appendix A Single site Penning trap

We consider for the moment a simplified setting, in which each ion is trapped in a symmetric static quadrupole potential $m\omega_{z}^{2}/2(z^{2}-(x^{2}+y^{2})/2)$ for ions of mass $m$ and charge $e$ embedded in a magnetic field of strength $B_{0}$ aligned along the $z$ axis. At a single site, the potential and magnetic field give rise to a Hamiltonian

[TABLE]

where $\omega_{1}=\sqrt{\omega_{c}^{2}-2\omega_{z}^{2}}$ and $\omega_{c}=eB_{0}/m$ is the bare cyclotron frequency. Writing the position and momentum operators in terms of creation and annihilation operators for the individual $x$, $y$, and $z$ degrees of freedom, defined as

[TABLE]

the Hamiltonian can be re-written as

[TABLE]

This can be separated into a sum of three independent harmonic oscillators using the transformation $\hat{a}_{\pm}=\frac{1}{\sqrt{2}}\left(\hat{a}_{x}\pm i\hat{a}_{y}\right)$ for the radial motion. We then obtain the expression given in equation 1 in the main manuscript.

Appendix B Normal modes: classical description

We consider a system of $N$ micro-Penning traps containing a single ion each (with charge $+e$) arranged arbitrarily in space. The Coulomb interaction between ions leads to a coupling between their motional states, resulting in $3N$ collective normal modes of motion.

Lagrangian Formulation

Let the quadrupole center $j$ and the position of the ion $j$ in the reference frame of the lab be defined by the coordinates $\mathbf{D}_{j}$ and $\mathbf{R}_{j}$ respectively. Then the local coordinates of the ion $j$ with respect to this quadrupole center are given by the vector $\bar{\mathbf{r}}_{j}=\mathbf{R}_{j}-\mathbf{D}_{j}$.

The trapping electrodes create a static quadrupole electric potential centered at each site $j$ and this potential can be written in terms of the local coordinates as $\phi_{j}=\sum_{\mu\nu}\phi_{j0}^{\mu\nu}\bar{r}_{j}^{\mu}\bar{r}_{j}^{\nu}$, where the indices $\mu$ and $\nu$ run over the Cartesian components, $x$, $y$ and $z$.

The electrostatic potential acting on the ion $j$ due to the Coulomb interaction with other ions is

[TABLE]

where $k_{e}=1/(4\pi\epsilon_{0})$ is the Coulomb constant.

The total electric potential, in the absence of any oscillating fields, is thus given by $\Phi_{j}=\phi_{j}+\kappa_{j}$.

A static homogeneous magnetic field $\mathbf{B}=B_{0}\sin\theta\cos\varphi\,\hat{\bf e}_{x}+B_{0}\sin\theta\sin\varphi\,\hat{\bf e}_{y}+B_{0}\cos\theta\,\hat{\bf e}_{z}$ creates the vector potential $\mathbf{A}_{j}$ at the site $j$. In the symmetric gauge, $\mathbf{A}_{j}=\frac{1}{2}(\mathbf{B}\times\mathbf{R}_{j}).$

In the laboratory frame of reference, the total Lagrangian of the system is then given by

[TABLE]

where $m_{j}$ is the mass of the $j$th ion.

The normal mode analysis begins by finding the equilibrium configuration of ions $\{\mathbf{R}_{j0}\}$, which is determined by the minimum of the total potential energy. By expanding the system Lagrangian about the equilibrium position of each ion in a Taylor series up to second order, we get a Lagrangian in terms of the generalized position vectors $\mathbf{r}_{j}=\mathbf{R}_{j}-\mathbf{R}_{j0}$ which specify the displacement of each ion from its equilibrium point. The second order term in the expansion effectively determines the normal mode dynamics of the system near the stable spatial configuration and is given by

[TABLE]

We proceed by putting together all the generalized position coordinates into a single $3N$-dimensional vector $q=\begin{bmatrix}x_{1}\,...\,x_{N}&y_{1}\,...\,y_{N}&z_{1}\,...\,z_{N}\end{bmatrix}^{T}$ and introducing the $3N\times 3N$ block matrices $M$, $W$, $V$ and $K$ constructed in terms of $N\times N$ sub-matrices as

[TABLE]

Here $\mathbb{I}_{N}$ and $\mathbb{O}_{N}$ represent the $N\times N$ identity and zero matrices respectively and the components of other sub-matrices are defined as

[TABLE]

[TABLE]

[TABLE]

where indices $j$ and $k$ run from $1$ to $N$ while again the indices $\mu$ and $\nu$ refer to the components $x$, $y$ and $z$.

The above definitions together with $\Phi=V+K$ allow us to write the effective phonon Lagrangian compactly as

[TABLE]

It should be clear that $M$ is a real diagonal matrix while $W$ is a real antisymmetric matrix. The matrix $V$ is traceless as a direct consequence of Laplace’s equation, while the matrix $K$ is traceless because the Coulomb forces being internal forces in the system of ions pairwise cancel each other and the total sum equates to zero. $V$ and $K$ are also both real and symmetric. As a result $\Phi=V+K$ is a real symmetric traceless matrix. These properties are useful in determining certain characteristics of the resulting normal mode eigenfrequencies and eigenvectors of the system.

Equations of Motion

Through the Euler-Lagrange equations,

[TABLE]

we can derive from the Lagrangian the equations of motion of our system. The two relevant derivatives are

[TABLE]

and

[TABLE]

which we can combine to get

[TABLE]

In vector form, we can then see that the equations of motion can be written as

[TABLE]

To find the normal modes of motion, we substitute the oscillating trial solution $q=q_{0}e^{-i\omega t}$ which yields the QEP

[TABLE]

that can be solved for complex eigenvectors $q_{0}$ and eigenvalues $\omega$, which in general can be complex. The set of eigenvalues $\{\omega_{\lambda}\}$ are the normal mode frequencies while the corresponding normalized eigenvectors $\{q_{\lambda}\}$ give us the normal mode coordinates.

The general solution can be written as

[TABLE]

where $\rho_{\lambda}$ are complex scalars. The motion of the ions in terms of the normal modes can then be retrieved as

[TABLE]

For real frequencies, the collective motion is bounded and hence all ions are confined.

Appendix C Normal modes: quantum mechanical description

From the Lagrangian of the system we can identify canonical conjugate variables to formulate our Hamiltonian. The generalized momentum corresponding to the generalized position $q_{j}$ is given by $p_{j}=\frac{\partial L}{\partial\dot{q}_{j}}$. The classical Hamiltonian of the system is then

[TABLE]

Quantizing the generalized coordinates as operators satisfying the standard commutation relations

[TABLE]

we can formulate the quantum mechanical Hamiltonian of the system as

[TABLE]

where $T=WM^{-1}W$ is a real symmetric matrix.

To diagonalize the Hamiltonian in the second quantized form $\hat{H}=\sum^{3N}_{\lambda=1}\hbar\omega_{\lambda}(\hat{a}^{\dagger}_{\lambda}\hat{a}_{\lambda}+\frac{1}{2})$, we form the phonon creation and annihilation operators, $\hat{a}^{\dagger}_{\lambda}$ and $\hat{a}_{\lambda}$, for the mode $\lambda$ as linear combinations of the generalized position and momentum operators,

[TABLE]

[TABLE]

where $\alpha_{\lambda k}$ and $\beta_{\lambda k}$ are complex numbers. For the commutation relation $[\hat{a}_{\lambda},\hat{a}^{\dagger}_{\lambda^{\prime}}]=\delta_{\lambda\lambda^{\prime}}$ to hold, the Hamiltonian must satisfy the commutation relation

[TABLE]

This commutator can be calculated by substituting $\hat{H}$ and $\hat{a}^{\dagger}_{\lambda}$ in terms of the canonical variables and comparing the coefficients of $\hat{p}_{l}$ and $\hat{q}_{l}$ in $[\hat{H},\hat{a}^{\dagger}_{\lambda}]=\hbar\omega_{\lambda}\sum^{3N}_{l=1}(\alpha_{\lambda l}\hat{p}_{l}+\beta_{\lambda l}\hat{q}_{l})$ yields the following set of coupled equations

[TABLE]

These can be written more succinctly in vector form as

[TABLE]

On eliminating $\beta_{\lambda}$ using $\beta_{\lambda}=i\omega_{\lambda}M\alpha_{\lambda}+\frac{1}{2}W\alpha_{\lambda}$, we then see that

[TABLE]

which is the same QEP encountered in the classical analysis in Appendix B. The QEP yields $6N$ eigenvectors and $6N$ eigenvalues. $3N$ eigenvectors will be used to form the creation operators while the other $3N$ eigenvectors to form the annihilation operators. We note that if the pair $(\nu_{\lambda},u_{\lambda})$ satisfies the QEP then the pair $(-\nu_{\lambda},u^{*}_{\lambda})$ also satisfies the QEP. Thus the total set of $6N$ eigenpairs

[TABLE]

for $\lambda$ running over $1$ to $6N$ can be divided into two equally sized subsets depending on the signs of the eigenvalues:

[TABLE]

The index $\lambda$ now runs from $1$ to $3N$ so that $\nu_{\lambda}$ is assumed to be positive from hereon.

Selecting the $3N$ eigenpairs which form the creation operators effectively means picking the sign of the eigenfrequency (and the corresponding eigenvector) for a given mode $\lambda$ in $\hat{H}=\sum^{3N}_{\lambda=1}\hbar\omega_{\lambda}(\hat{a}^{\dagger}_{\lambda}\hat{a}_{\lambda}+\frac{1}{2})$ and involves fixing the normalization of the eigenvectors $\alpha_{\lambda}$ so that $\left[\hat{a}_{\lambda},\hat{a}^{\dagger}_{\lambda}\right]=1$. Explicitly,

[TABLE]

Substituting $\alpha_{\lambda}=c_{\lambda}\gamma_{\lambda}$, where $\gamma_{\lambda}$ is normalized to one and $c_{\lambda}$ is a complex scalar,

[TABLE]

which for the condition $\left[\hat{a}_{\lambda},\hat{a}^{\dagger}_{\lambda}\right]=1$ yields

[TABLE]

Since $|c_{\lambda}|^{2}$ is non-negative, $\left[\hat{a}_{\lambda},\hat{a}^{\dagger}_{\lambda}\right]=1$ only when the quantity $\omega_{\lambda}/(\omega^{2}_{\lambda}\gamma^{H}_{\lambda}M\gamma_{\lambda}+\gamma^{H}_{\lambda}\Phi\gamma_{\lambda})$ is positive. This expression helps us pick out the $3N$ eigenpairs to form the creation operators

[TABLE]

where

[TABLE]

can without loss of generality be chosen as real. Further, inverting the expressions for the creation and annihilation operators yields the second quantized form of the position and momentum operators,

[TABLE]

and

[TABLE]

Appendix D Trap imperfections

In a real experimental setup, the trapping potential may not be of the idealized form expected from the optimization of the electrode structures, while the magnetic field could be misaligned with the confining direction of the potential. As long as the imperfections in the electric potential are harmonic and the magnetic field is homogeneous over the entire system, we could employ the general discussion in Appendix B in order to study the normal modes of the imperfect system.

Linearization of the QEP 40 in the first-companion form yields the GEP

[TABLE]

Inversion of the matrix on the r.h.s. of the above equation leads to a further reduction to the SEP

[TABLE]

with $6N$-dimensional eigenvectors $v=\begin{bmatrix}q_{0}&\omega q_{0}\end{bmatrix}^{T}$ and $6N$ eigenvalues $\omega$ belonging to the $6N\times 6N$ matrix

[TABLE]

For the sake of brevity, we define the matrices $W^{\prime}=M^{-1}W$ and $\Phi^{\prime}=M^{-1}\Phi$, so that we have

[TABLE]

and

[TABLE]

Since $A^{2}v=\omega^{2}v$ and the sum of eigenvalues of a matrix is equal to its trace,

[TABLE]

where we use the fact that $\Phi^{\prime}$ is traceless. The stability of the system can as usual be determined by checking if all eigenvalues are real. Noting that the frequencies come in pairs of positive and negative values in the stable regime we can express the above sum in terms of the $3N$ positive frequencies,

[TABLE]

This trace can be explicitly calculated as

[TABLE]

where $\omega_{c,j}=eB_{0}/m_{j}$ is the true cyclotron frequency of the $j$th ion, thus allowing us to express the sum as

[TABLE]

Note that the sum is independent of the trapping potential. For a typical experiment with ions of the same species and no impurity defects, $m_{j}=m$, and the above sum further simplifies in terms of the common true cyclotron frequency $\omega_{c}=eB_{0}/m$ to

[TABLE]

Equation 70 can be treated as a non-trivial generalization of the well known Brown-Gabrielse invariance theorem for a single ion in an imperfect Penning trap,

[TABLE]

One additional result can be derived from equation 62 by using the property that the product of the eigenvalues of a matrix is equal to its determinant so that

[TABLE]

An interchange of $3N$ columns in the matrix on the r.h.s. allows us to write the product in terms of the $3N$ positive frequencies as

[TABLE]

or

[TABLE]

Finally, we arrive at

[TABLE]

which for the case of ions having identical masses can be more conveniently expressed as

[TABLE]

This result tells us that the product of eigenvalues is independent of the magnetic field and depends only on the curvature tensor of the total electric potential.

Appendix E Spin spin coupling

The derivation in this Appendix follows closely the methodology from ref. Wang et al. (2013). The ODF leads to the interaction term

[TABLE]

In the Lamb-Dicke regime, we can expand this expression in terms of the equilibrium positions and deviations from them as

[TABLE]

Then the effective spin Hamiltonian is given by

[TABLE]

which uses the definitions

[TABLE]

[TABLE]

and

[TABLE]

We can express the excursions from equilibrium in terms of the harmonic oscillator creation and annihilation operators, giving

[TABLE]

In the interaction picture with respect to the oscillator mode Hamiltonian $\hat{H}_{\text{PH}}$, we then find that

[TABLE]

where we have defined the functions

[TABLE]

From equations 79, 80 and 81, and making the rotating wave approximation with respect to the oscillator frequencies, we then find that the static part of the effective spin Hamiltonian can be written in the form of an Ising-like spin Hamiltonian

[TABLE]

with the coupling terms given by

[TABLE]

Here $\gamma_{\lambda}$ is the normalized normal mode eigenvector corresponding to the frequency $\omega_{\lambda}$ and the indices $\nu,\nu^{\prime}$ run over $x,y,z$.

Appendix F Quantum gates

The theory for calculations of quantum gates is similar to that detailed above, with the difference that for fidelity calculations we have to consider the effects of residual spin-motional entanglement and geometric phases from motional state components displaced in phase space. For a gate operated using an oscillating force with state-dependence in the $z$-basis, a common method for measuring the fidelity is to sandwich the two-ion gate in a Ramsey experiment performed simultaneously on both ions. This produces the maximally entangled state $\left|{\psi_{\rm B}}\right>=\left(\left|{00}\right>-i\left|{11}\right>\right)/\sqrt{2}$. Thus the force pulse acts on a superposition state $\left(\left|{00}\right>+\left|{01}\right>+\left|{10}\right>+\left|{11}\right>\right)/2\equiv\sum_{E}\left|{E}\right>/2$ with $E\in\{00,01,10,11\}$. If the motional state of all modes is prepared in the ground state, the initial state of internal and motional states is

[TABLE]

For two ions with index $j=1$ and $j=2$ which are at the same phase of the optical dipole force (here assumed to be zero), the interaction picture Hamiltonian $\hat{V}_{I}(t)$ becomes

[TABLE]

where $\rm h.c.$ indicates the Hermitian conjugate. Since the gate acts on the four eigenstates of $\sigma^{z}_{1}+\hat{\sigma}^{z}_{2}$, the action of the Hamiltonian acts in terms of these four eigenstates as

[TABLE]

where $\hat{D}_{\lambda}$ denotes a displacement operator on the mode $\lambda$, with the displacement amplitude

[TABLE]

with the function

[TABLE]

and $p_{\lambda,E}=\sum_{\nu}k_{R}^{\nu}\left(\pm\alpha_{\lambda 1\nu}+\pm\alpha_{\lambda 2\nu}\right)$ with plus signs for internal state [math] and minus signs for internal state $1$ of each ion. The phases are given by

[TABLE]

with

[TABLE]

Once the state after the gate has been found, it can be used to form a density matrix from which all relevant quantities can be obtained. For the fidelity with the maximally entangled state, we take $F=\left<{\psi_{B}}\right|\rho(t)\left|{\psi_{B}}\right>$. We trace out the motional displacements using the overlap between two coherent states

[TABLE]