Quantum simulation of the hexagonal Kitaev model with trapped ions

We present a detailed study of quantum simulations of coupled spin systems in surface-electrode ion-trap arrays, and illustrate our findings with a proposed implementation of the hexagonal Kitaev model [A. Kitaev, Annals of Physics 321,2 (2006)]. The effective (pseudo)spin interactions making up such quantum simulators are found to be proportional to the dipole-dipole interaction between the trapped ions, and are mediated by motion which can be driven by state-dependent forces. The precise forms of the trapping potentials and the interactions are derived in the presence of a surface electrode and a cover electrode. These results are the starting point to derive an optimized surface-electrode geometry for trapping ions in the desired honeycomb lattice of Kitaev's model, where we design the dipole-dipole interactions in a way that allows for coupling all three bond types of the model simultaneously, without the need for time discretization. Finally we propose a simple wire structure that can be incorporated in a microfabricated chip to generate localized state-dependent forces which drive the couplings prescribed by this particular model; such a wire structure should be adaptable to many other situations.


Introduction
Ion trap systems have proven to perform well for implementing the basic elements of traditional quantum computing, where evolution is described in terms of discrete gate operations, which can be implemented step by step as intermediate states are irrelevant. This is in contrast to quantum simulations, where the goal is to simulate the continuous evolution of a given Hamiltonian. While the initial proposal for quantum computing with trapped ions relied on a number of sequential steps to mediate effective qubit interactions [1], other approaches [2][3][4][5][6][7][8][9][10] achieve interaction between the internal states of the ions via constant Hamiltonians and therefore allow the development of quantum simulators based on trapped ions [7,8,[11][12][13][14][15]. In such simulators, interactions between trapped ions are dominated by the Coulomb potential. For this interaction to affect internal states (i.e., the qubits or pseudo-spins representing the effective quantum system to be simulated), state-dependent forces must be applied to some or all of the trapped ions. State-dependent forces can be achieved through optical ac Stark shifts [5,7,[15][16][17][18][19], static magnetic-field gradients in combination with homogeneous radio-frequency (rf) fields [4,10,14,20], or with rf field gradients [21,22]. While in most cases the Coulomb interaction is considered between ions in a self-assembled single chain or crystal, coupling of independently trapped ions has recently been demonstrated [23,24].
For quantum simulations with ions in microtraps, we must take into account how the presence of the electrodes modifies the Coulomb interaction. While in many systems this effect is negligible (for example, in the surface-electrode setup of [23] the Coulomb coupling was found to be enhanced by only 1.8 %, in agreement with our more general results in section 2.3), the general theory developed here for a lattice of surface-electrode (SE) microtraps shows that significant modifications to free-space couplings are possible. Far from being an inconvenience, these modified interactions can be used to design quantum simulations with specific short-range effective pseudo-spin interactions, which we illustrate with the hexagonal Kitaev model as a concrete example.
The remainder of this paper is organized as follows. In section 2 we present a Green's function approach to solving electrostatic problems as they occur for surfaceelectrode ion traps in the presence of a cover electrode. In section 3 we derive general expressions for spin-spin couplings in two-dimensional microtrap arrays, applicable, for example, to electric coupling to light fields or magnetic coupling to microwave nearfield gradients. In section 4 we combine all these methods to show how the hexagonal Kitaev model [25] can be implemented with an array of trapped ions on an optimized surface-electrode chip, including a dedicated wire structure that could be integrated in the chip to simultaneously mediate the couplings along three distinct bonds by use of magnetic-field gradients. Finally, Appendix A gives a summary of the used coordinate systems.

Electrostatics in the presence of conducting planes
The electrostatic interaction between charged particles close to conducting surfaces can be strongly modified by the presence of the conductors [27]. In the idealized geometry of a perfectly conducting grounded electrode plane at z = 0, the total electrostatic energy of a set of charges Q i located at positions r i in the half-space z i > 0 (see figure 1) is [26] The Coulomb interaction term in (1) is expressed in terms of the Dirichlet Green's function G ∞ , which can be found from the free-space Green's function G (0) (r, r ) = 1/ r − r by the method of images (see figure 1),  Figure 1. Coulomb interaction between two charges Q and Q (full red circles) in the presence of a grounded plane. The image charges (empty red circles) are located below the grounded plane and carry opposite charge. Interactions between the charges (full blue arrow) contribute to (1) in full, while interactions between charges and image charges (dashed blue arrows) contribute with a prefactor of 1/2 [26].
where ρ = (x − x ) 2 + (y − y ) 2 is the horizontal distance between the charges. In the following we review the effects of a grounded cover plane, i.e., a second, parallel conducting plane covering the electrode plane at height z = H (see figure 2a). In the initial proposal [28] and demonstration [29] of surface-electrode rf traps, the conducting surface nearest to the trap electrodes was theoretically at infinity but in practice a part of the surrounding apparatus. It has been suggested that adding a cover plane in the form of a dc-biased mesh above the electrodes could improve trap depth [30]. In addition to possible benefits of providing bias field and shielding, the cover plane could have more practical advantages, namely shielding the trapping region from fields due to quasi-static charges on insulators in the vacuum chamber, and establishing a more well-defined boundary condition. Further, if the cover plane is modified to carry rf and dc electrodes of arbitrary shape in the same way as the electrode plane, the presented formulas can be used directly to calculate the combined electric fields generated in this "sandwich trap" geometry (however, if optical access to such a trap geometry is achieved with holes and/or fiber optics in the electrode planes [31], the present full-plane treatment must be adapted [32]).
Below, we first modify the Green's function (2) to include the cover plane and illustrate that a cover plane at height H leads to exponential shielding on a lateral length scale of H (section 2.1), then consider its effects on the electric potential generated by surface electrodes (section 2.2) and on effective dipole-dipole interactions between vibrating trapped ions (section 2.3).

the shielding effect of the cover plane
When a grounded conducting cover plane at height z = H is added to the setup of figure 1, the Coulomb interaction (1) of charges located between these two planes is modified to (a) Sketch of a surface-electrode trap with a grounded cover plane positioned at a height H above the electrode plane. The red ring electrodes are at rf potential, while all grey areas are grounded. For static interactions or interactions varying slowly compared to the rf period, only the time-averaged potential contributes, so for our purposes the situation is equivalent to two completely grounded planes. (b) The interaction energy (3) between two point charges at same height h over the electrode plane, as a function of the charge separation ρ, in the presence of a cover plane at height H = 100h. The red and blue parts of the solid curves are computed by (4a) and (4b), respectively, while the dashed lines illustrate the approximate behavior given by (5).
Both the scaled self-potential e H (z) and the Dirichlet Green's function G H corresponding to the cover plane geometry with infinite conducting electrode planes at z = 0 and z = H can be found by summing over an infinite sequence of mirror planes; and in the absence of a cover plane (H → ∞) they reduce to (1) and (2). The scaled self-potential is in terms of Euler's constant γ = 0.577216 . . . and the digamma function ψ(a) = Γ (a)/Γ(a). The Dirichlet Green's function is where K 0 is the modified Bessel function of the second kind. The second form (4b) is obtained by solving the Laplace equation in cylindrical coordinates [27]. Both forms converge for all parameters (ρ, z, z ); but while (4a) converges faster when r − r H, (4b) is more suitable if r − r H, in particular for ρ H, as discussed below. The Coulomb interaction energy G H (ρ, z, z )QQ /(4π 0 ) between two charged particles in a SE trap depends on the horizontal separation ρ, as illustrated in figure 2b. To illustrate this interaction energy we take the particles to be at the same height h above the electrode plane, and the cover plane height H to be much larger than h. When ρ H, we expect the cover plane to be irrelevant, so that the interaction is described by a single image charge: it falls of as ρ −1 while ρ < h (where the electrode plane is irrelevant) and as ρ −3 thereafter, as described by (2). When ρ H the cover plane becomes important and the asymptotically dominant form is the first term of the resummation (4b), so that the presence of the cover plane leads to an exponential shielding at the length scale of the cover plane height, as illustrated in Fig. 2. Summarizing, 2.2. the potential due to the surface electrodes The contribution to the total potential from the structured electrodes in the z = 0 plane can be computed as an integral over the electrode plane: where we have introduced a "surface Green's function" G (S) In the absence of a cover plane, the surface Green's function was found to be [32] with the geometric interpretation that the potential at r due to an electrode at potential Φ 0 is Φ 0 /2π times the solid angle spanned by the electrode as seen from r [33]. Alternatively, the electric field at r is proportional to the magnetic field that would be observed if a current were flowing along the edge of the electrode [33,34]. For electrode configurations that are translationally invariant in the x-direction, the system can be described by conformally mapping the upper-half yz-plane (z > 0) to a disc [35]. Analogous to (4a) and (4b), we have for the general case including a cover plane, where ζ is the Riemann zeta function, P j are Legendre polynomials, and s = r − r = ρ 2 + z 2 . Forms (8a) and (8b) converge for all (ρ, z); but while (8a) converges faster when s H, (8b) is more suitable if s H. Form (8c) is restricted to s < 2H and is most useful for s H. Similar to (5) we find the approximate behaviors where ψ (a) = Γ (a)/Γ(a) − ψ 2 (a) is the first derivative of the digamma function. We conclude that the influence of any surface electrode is exponentially damped at distances larger than H, which is advantageous for the experimental construction of quasi-infinite surface microtrap lattices in that it reduces the influence of the inevitable electrode boundary: at any point further than H away from the edge of the electrode and cover plane, the trap will look as if the electrode were infinitely large.
Since the surface Green's function only depends on the x and y coordinates through r − r , (6) is a folding integral (convolution) [32] and can be rewritten as a product of the Fourier-transformed quantities,Φ(k x , k y , z) =G and a similar expression for the Fourier-transformed Green's function. The latter is cylindrically symmetric (k = k 2 x + k 2 y ), and allows a rather intuitive interpretation. All solutions of the Laplace equation with horizontal wavevector {k x , k y } are of the form e i(kxx+kyy) (α + e +kz + α − e −kz ); the Green's function (11) gives the unique solution which has unit amplitude on the electrode plane [G . Therefore (11) gives the unique extension of a unit-amplitude potential plane wave from the z = 0 plane into the z > 0 half-space which satisfies the boundary condition of vanishing amplitude on the cover plane. The fact that the momentum-space representation of the surface Green's function (11) can be written without infinite sums greatly simplifies the description of infinite lattices of surface-electrode microtraps [36].

dipole-dipole interactions between trapped ions
Trapped-ion quantum simulators couple internal degrees of freedom of the ions (typically hyperfine states or metastable D-states) through a state-dependent coupling to shared vibrational degrees of freedom [1,6,8,14,20,21] (see section 3). A crucial ingredient of these couplings is the precise nature of the Coulomb interactions between the ions. Here we address the details of this latter point, since it will determine how to construct a quantum simulator of a desired system, as exemplified in section 4.
We consider the regime of "stiff" ion trapping [6], where the Coulomb interaction is relatively small compared to the trapping potential, and we can interpret the normalmode dynamics of the ion crystal as that of a set of local harmonic oscillators that are weakly coupled. The ion trapping potential defines a set of local eigenmodes for the i th ion corresponding to vibration in three orthogonal directions m µ i (with m µ i = 1 for µ = 1, 2, 3) around an equilibrium position R 0i . In what follows we use these directions to parametrize the position of the i th ion as The total Coulomb energy of a set of N charges is given in (3), and the leading-order term that couples the motion of the ions is Since we are mainly interested in near(est)-neighbor interactions, we evaluate this expression in terms of the infinite sum over image charge pairs (4a), rather than the resummed form (4b): where the explicit dipole-dipole coupling is given by the expression without a cover plane, in terms of n = (r − r )/ r − r ,n = (r −r )/ r −r , and the mirrored quantities r = r − 2(r ·ẑ)ẑ andm = m − 2(m ·ẑ)ẑ. The first term of (15) is the well-known dipole-dipole interaction, while the second term is the correction due to image charges in the electrode. In order to illustrate the behavior of the dipolar interaction (15) in close proximity of a conducting electrode plane, we again consider two ions located at equal height h above the electrode plane, spaced by a distance ρ along the x axis, and in the absence of a cover plane. If we assume that both ions vibrate along axes m = m that are parallel to the lab-frame coordinate axes, then we find that the presence of the electrode plane can either increase or decrease the dipolar coupling strength: Thus we see that by choosing the directions of vibration m in particular ways we can use the presence of the electrode plane to make the dipolar interactions fall off with the fifth power of distance instead of with the third power (which is the case in the absence of any conducting planes), as long as the ion oscillation frequency is low enough to avoid the effects of retardation and dissipation. The relevant length scale that determines whether or not the electrode plane has a strong influence on the dipole-dipole coupling is ρ ∼ 2h, similar to figure 2; for even farther separations (ρ > H) we find exponentially damped dipole-dipole couplings due to the shielding effect of the cover plane (see section 2.1). These rapid dampings can be used to construct lattice simulation models with nearly local interactions, which is a desirable feature since many spin models from condensedmatter physics are formulated in terms of such local (e.g., nearest-neighbor) couplings.

Spin-spin interactions between trapped ions
This section derives how state-dependent forces can induce pseudo-spin interactions between neighboring ions through the Coulomb potential. While this effect is wellknown in principle [1], we show how these effective interactions are constructed in a lattice of ions without the need for time-slicing ("Trotterization" [37]). Further we show that, to lowest order, the effective interaction strengths are proportional to the realspace Coulomb coupling strengths, an observation that greatly simplifies the design of lattice-based quantum simulators (see section 4).

normal modes of vibration
For small oscillation amplitudes r µ i the coupled harmonic motion of N ions in a lattice can be described by considering the local trapping potential curvatures (the second derivatives of the ion trapping pseudopotential with respect to position) around the equilibrium positions R 0i and including the Coulomb couplings between ions in separate wells to second order [38]. If we assume that all excursions r i − R 0i are already written in terms of the "bare" eigenmodes of the isolated local trapping potentials (with "bare" frequenciesω iµ ), as in (12), then the potential energy of the ions of mass M is where γ µν ii = 0 and γ µν i =j = (13). This quadratic potential energy can be diagonalized using coefficients O iµm such that the real-space displacements can be written as and V = 3N m=1 1 2 M ω 2 m q 2 m in terms of the lattice normal-mode amplitudes q m and their frequencies ω m . We quantize these normal modes through q m →q m = q 0m (â m +â † m ) with q 0m = /(2M ω m ) and the usual commutation relations [â m ,â † m ] = δ mm . We will work in the "stiff" lattice limit, where we assume that the "bare" trap frequencies ω iµ ≡ω µ ∀i are equal for all ions ‡ and the Coulomb couplings between ions do not significantly mix local modes with different values µ. This is the case when (i) the bare trap frequencies are sufficiently far apart: |γ µν ij | |ω 2 µ −ω 2 ν | ∀i, j, µ = ν, and (ii) the vibrational bands are sufficiently narrow: |γ µµ ij | min ν |ω 2 µ −ω 2 ν | ∀i, j, µ. In this limit, the normal modes of the lattice separate into three disjoint sets indexed by µ ∈ {1, 2, 3}, each containing N normal modes with frequencies close to the correspondingω µ . ‡ This equality can be relaxed to the condition that for the modes that are used for inducing spin-spin couplings, the spreads of the bare frequencies are much smaller than the dominant Coulomb couplings.

state-dependent forces
In addition to the Coulomb couplings, which are always on and define the coupled vibrational eigenmodes of the trapped ions, we can experimentally introduce fields that couple to internal states of the ions. Examples of such interactions are electric or magnetic dipole couplings, Raman couplings, or electric quadrupole couplings to laser or microwave fields. In the following general treatment we assume that there is a coupling between (a) classical external field(s) effectively oscillating with angular frequency ω I , and two internal states of each ion, forming an effective two-level (spin-1/2 or pseudospin-1/2) system. Irrespective of the type of induced coupling, the coupling operator of the i th ion in its (pseudo)spin-1/2 subspace can be expressed as a linear combination of the identity operatorσ where we have performed a first-order expansion in the ion positions assuming small oscillation amplitudes. Any type of spin-1/2 coupling that is used with trapped ions (including effective couplings to pseudo-spin degrees of freedom) can be brought into this form, where terms with non-vanishing prefactors c (i) and s (i) are referred to as "carrier" and "sideband" terms, respectively. The phases can absorb differences in the details of driving fields: while stationary fields in general have φ s , travelling waves (e.g., light fields) are characterized by φ As an example, the coupling of physical spins to a magnetic field is found by expanding their magnetic dipole operators asμ ZẐ ) in terms of the Bohr magneton µ B and the g-factors g (i) ; for small ion excursions the coupling Hamiltonian to the magnetic field H I = − iμ (i) · B(r i ) cos(ω I t + φ) can thus be expressed in the form of (19) with for = X, and similarly for = Y, Z and φ We stress, however, that the above form of the magnetic dipole operator does not apply to pseudo-spins for their effective interactions with external fields; see section 4.4 for an example involving pseudo-spins.

effective spin-spin interactions
Inserting the lattice normal-mode expansion (18) and (12) into (19), we can write the interaction Hamiltonian as where we have dropped the approximation symbol and introduced Ω im = It is common to transform into the interaction picture to assess the dynamics induced by such an interaction Hamiltonian. In this transformation, the field-free Z leads to a time dependence of the operators in (21): The terms involvingσ Y can lead either to spin flips without affecting the motion ("carrier"-transitions, mediated by c  Y and resonant around ω ↑↓ ±ω µ ); the latter will dominate if they are not driven too strongly and |ω I − ω ↑↓ ±ω µ | |ω I − ω ↑↓ | for one of the signs in ±. Here we concentrate instead on a drive with frequency |ω I −ω 3 | ω 3 ω ↑↓ , close to one of the three bare eigenfrequencies of the uncoupled ion sites (we have chosen µ = 3 without restricting generality). In this case we can neglect the terms in s after a second rotating-wave approximation, with the detunings δ m ≡ ω m − ω I . Equation (23) can be exactly integrated via a Magnus expansion [39,40] to yield the unitary evolution operator s . The first exponent describes a set of time-dependent coherent displacements to all normal modes that can entangle the motion with the internal pseudo-spin states. The second exponent constitutes a phase that depends on pairs of spins and can be interpreted as a spin-spin interaction. For a faithful simulation of interacting spins it is desirable that (A) the first term should be as close as possible to the identity operator, in order to avoid populating vibrational excitations, and (B) the second term should provide sizable phases for desired inter-ion couplings with i = j, as these represent the spin-spin interactions. It can be shown from the expression above or by the use of a canonical transformation [6] that (A) can be approximately met as long as |Ω im | |δ m | for all (i, m, ). The above restrictions do not limit the time scale for simulations, as long as one assumes that sufficiently strong couplings can be induced by lasers or microwave field gradients. However, the energy scale of nearest-neighbor Coulomb interactions also plays an important role in determining simulation time scales, but this dependence is hidden in the normal-mode coefficients O iµm of (18). To illustrate this point, we assume that Ω im0 = 0 for all normal modes m and sites i [see (20) for an example]; but what we show below also holds for more general cases. Assuming negligible displacements [see point (A) above] the unitary evolution operator thus simplifies tô In the "stiff" lattice limit (see page 8) we choose the drive frequency ω I close to one of the bare frequencies, sayω 3 , such that the detunings δ m will be much smaller for normal modes in this set than for the other normal modes; consequently, the sum over modes m in (25) can be restricted to an "active" set of N modes clustered aroundω 3 . If we further choose the drive frequency such thatδ 3 =ω 3 − ω I is much larger than the spread of the normal mode frequencies in the active group, the series expansion together with the relations in the "active" group of normal modes simplifies the unitary evolution (25) tô where we have further approximated q 0m ≈q 03 = /(2Mω 3 ). The first term in (28) is a global phase; it is the second term that mediates an effective spin-spin coupling on the lattice of ions. It can be interpreted as the evolution under an effective spin-spin coupling Hamiltonian due to the driving of local mode µ = 3, where the effective spin-spin coupling coefficients are We conclude that to lowest order the strength of the spin-spin coupling between two ions is determined by the geometric overlaps (m 3 i · s (i) Z ) of the "active" local modes of vibration with the direction of the state-dependent force, as well as by the real-space Coulomb coupling strength γ 33 ij between the ions moving along these local modes [see (13) and (17)]. The relative phases φ ij s of the driving forces can be used to modulate the coupling strengths. These observations are used in section 4 to construct a quantum simulator on a lattice of trapped ions. Equations (29) and (30) faithfully describe the evolution of the system under the following conditions: • The vibrational band structure of the trapped ions must consist of clearly distinct bands which can be addressed individually (see page 8 for the conditions for "stiff" trapping).
• The state-dependent force must be driven at a small detuning δ from one of these bands: δ must be large enough such that (26) is valid for this band, but small enough such that the contributions to (30) from other bands are negligible.
• The amplitude of the state-dependent force must be small enough such that it does not significantly excite ion vibrations. The condition |Ω im | |δ m | given above implies that the state-dependent forces must be weaker than the force scalê F ∼ |δ|/q 0 for the addressed band.
The above arguments can be made in a very similar fashion forσ interactions by considering the interaction Hamiltonian (21) with ω I ≈ ω ↑↓ ±ω µ in the appropriate basis |± = (|↑ ± e iχ |↓ )/ √ 2, where the considered spin-spin interaction is diagonal. The only slight complication can arise from carrier terms proportional to c (i) X and c (i) Y that are detuned by roughly the motional eigenfrequencies. For detunings from the sidebands on the order of the dipole-dipole interactions and correspondingly small drive strengths, however, these carrier terms can be safely neglected.
To summarize this section, we have considered general effective spin-spin interactions in the limit of "stiff" ion trapping. We have shown that even in a lattice, the spin-spin coupling strength of any two ions depends on the dipolar Coulomb coupling between these two ions. To avoid appreciable entanglement between (pseudo)spins and ion motion, the detunings of driving fields need to be larger than the couplings they induce. This latter finding agrees with other work on simulation with trapped ions [8]. At the same time, the detunings cannot be much larger than the couplings between nearest neighbors, which determine the finer structure of the normal-mode spectrum around the frequencies of the uncoupled ("bare") motion of an ion tightly bound in one of the trapping wells (see section 4 for a concrete example). These requirements impose stringent bounds on the time scales necessary to perform simulations. For example, in [23] two ions at a distance of 40 µm exhibited an exchange splitting of approximately 3 kHz, barely sufficient to demonstrate a few energy exchanges before ion heating profoundly altered the motion. Simulations that need to progress adiabatically with respect to this exchange period will therefore be experimentally challenging and may require reducing anomalous heating below what was measured in [23].

Kitaev model
As an example of how to use the results of sections 2 and 3 in the design of a quantum simulator, we construct an implementation of the hexagonal Kitaev model [25] with microtrapped ions. In its ideal form, this exactly solvable two-dimensional spin model has a topologically ordered ground state with anyonic excitations, which makes it extraordinarily interesting for study in a quantum simulator with individual access to the constituent degrees of freedom.

model and implementation
The hexagonal Kitaev model [25] has the Hamiltonian defined on a honeycomb lattice of spin-1/2 particles, where the lab-frame bond vectors refer to Fig. 3: In this way the Hamiltonian (31) associates each real-space bond direction (32) with a spin quantization direction; however, it is important to keep in mind that the bond directions and the associated spin quantization directions are not a priori related (see Appendix A). The ions are located on two sublattices L • and L • , as shown in figure 3. Neighboring ions are a distance d apart.
The form of (31) is exactly that of (29) summed over three concurrent driving force fields. As these driving fields will be at very different frequencies, they can be applied simultaneously in order to drive the full Hamiltonian (31). What is therefore needed in order to implement the Kitaev model is a set of "bare" vibrational directions of the ions such that the couplings γ µν ij , and therefore the effective spin-spin couplings (30), match the particular geometry of the three terms in (31).
We choose the ion trapping height to be half of the inter-ion distance, h = d/2, and the orthonormal principal axes of vibration for ions on the two sublattices as  This particular choice of axes of vibration has the property that dipole-dipole couplings of the sort of γ µµ ij (i.e., coupling the m µ i vibration of the ion at R 0i with the m µ j vibration of the ion at R 0j ) are strongly dominated by the nearest-neighbor couplings required by the Kitaev model (31), shown in figure 3 for µ = X. These couplings can be calculated from (15) (in the absence of a cover plane); the resulting nearest-neighbor terms in the dipole-dipole coupling part of the Coulomb potential (17) are In addition, there are dipole-dipole couplings to neighbors that are further away and that turn out to be larger than the off-diagonal terms (µ = ν) in (34). The vibrational normal-mode band structure due to all of these dipole-dipole couplings is shown in figure 4, with the effective density of states shown in figure 5. It consists of two bands, in which neighboring ions oscillate in-phase (upper band) and out-of-phase (lower band), and whose small frequency spread is indicative of the dominance of the nearest-neighbor coupling over all other couplings. Many dipole-dipole couplings of the sort of r µ i r ν j with µ = ν are nonzero in this configuration; however they do not lead to effective spin-spin couplings if the underlying trap frequencies along the directions m µ i and m ν j are strongly off-resonant (see sections 3.3 and 4.2). Thus neglecting any µ = ν couplings, the effective spin-spin  ω P(ω) Hamiltonian that is constructed from µ = X Coulomb interactions is approximately where the first sum contains couplings between the different lattices while the second sum contains couplings within the lattices; numerical prefactors for the perturbing terms are used for brevity, as in figure 3. Here we have assumed for simplicity that all ions are simultaneously pushed by the same state-dependent force with equal phase. The Hamiltonians H Y and H Z are found from (36) through rotations by ±120 • , and the total effective spin Hamiltonian is where J X , J Y , and J Z are effective coupling constants containing the diagonal coupling strength, the physical prefactors, as well as the mechanisms used for achieving these effective spin-spin couplings (see section 4.4). The topic of whether or not this slightly perturbed Hamiltonian (36) exhibits the same interesting topological phases as the ideal Hamiltonian (31), at zero or finite [41,42] temperature, is beyond the scope of this article. We mention, however, that if the perturbative terms of (36) will be deemed too strong, they can be reduced further by driving the different wires with different relative phases or amplitudes [see (30)]. The presented configuration of trapping height and vibrational axes nearly maximizes the desired dipole-dipole couplings at the same time as it nearly mimimizes all undesired couplings. By numerical optimization we can identify a configuration that performs a few percent better than (33), but we have not been able to obtain an analytical description of this configuration.

surface-electrode trap design
To have maximally incommensurate vibrational frequencies along the normal mode axes (33) we choose them in the golden ratio ω X : ω Y : ω Z = φ −1 : 1 : φ with φ = (1 + √ 5)/2. We use the algorithm of Ref. [36] to find an rf surface-electrode pattern that will generate an infinite honeycomb lattice of exactly such microtraps, with the following constraints: • The unit cell of the electrode pattern is defined by the vectors a = d{ √ 3, 0, 0} and b = d{ √ 3/2, 3/2, 0}.
• The gradient of the rf electric potential generated by the surface electrodes must vanish at the ion positions in order to have minima of the rf pseudo-potential. The resulting electrode pattern is shown in the left panel of figure 6. It generates microtraps at the desired positions with dimensionless curvatures [36] κ = 0.080 and no spurious additional microtraps. This is to be compared with a simple out-of-plane quadrupole honeycomb lattice geometry (κ = 0.102) as in Ref. [36] (see figure 6, right panel), which can potentially be deformed during the experiment via dc electrode potentials into satisfying the above constraints. Such dc electrodes might be necessary in any experimental implementation in order to null micro-motion of the ions [43] induced by manufacturing inaccuracies and stray charges.

trap depth and trap loading
The depth of the microtrap lattices generated by the electrodes shown in figure 6 are rather shallow. For the honeycomb lattice, which is more easily analyzed due to its p6m symmetry, figure 7  adiabatic trap depth is likely further reduced by the breakdown of the pseudo-potential approximation near the trap barrier. In order to reliably load these microtraps we can make use of the cover plane (see section 2.1): applying a positive dc bias potential to the dc electrodes and the same potential to the cover plane adds an electrostatic potential that pushes the ions towards the rf electrode and increases the trapping well depth (see figure 7). Since this dc potential is equivalent to applying a negative dc bias potential to the rf electrode, its dc electric field at the ion trap sites vanishes (by construction of the rf electrode shape) and it thus does not induce micro-motion [43]. We find that applying a small dc bias potential of at least 0.125E pp /Q ≈ 0.6 V is sufficient to make the desired lattice of microtraps the only minima of the total potential (dashed line in figure 7). By applying a stronger bias voltage, the resulting total potential (dotted line in figure 7) is deep enough to trap ions produced by photoionization directly from a hot atomic beam. This bias will simultaneously cause the traps to be shallower in the xy plane.

wires for magnetic interaction
As described in section 3, effective spin-spin interactions between ions require internalstate-dependent forces to be applied to the ions. In the present model we propose to embed parallel wires below the electrode plane, which generate local magnetic field gradients at the positions of the ions [see (20)]. A relatively simple periodic grid of two different types of wires, indicated in red and blue in figure 3, suffices to implement the spin-spin interactions along all three bond types of the Kitaev model. As explained in section 3, one can induce pairwiseσ Y interactions by currents at frequencies that are near-resonant to ω ↑↓ ± ω X and ω ↑↓ ± ω Y (Mølmer-Sørensen type To avoid sizable entanglement of the (pseudo)spins with the ion motion, we need to fulfill The hexagonal Kitaev model features interesting gapped phases with anyonic excitations of the ground state for example if |J x | = |J y | < |J z |/2, in which case the coupling constant of the resulting effective Hamiltonian is J eff = J 2 x J 2 y /(16|J z | 3 ) < |J z |/256 {see [25] for discussions of these phases and their emergence from (31)}.
While the geometric prefactor of (40) depends on the details of the model, its functional dependences are expected to remain the same for a broad class of wire-driven coupled pseudo-spin models, in particular also for J X and J Y of the same system. The exact form of the interactions along the X and Y bonds depends on the transition dipole matrix elements µ d = ↑|μ|↓ which can have components along all spatial directions. § The component along the quantization axisẐ is relevant for π transitions (where thê Z component of the total angular momentum F of the ion does not change, ∆m F = 0) while perpendicular components can be used for σ ± transitions (during which theẐ component of F changes by ∆m F = ±1). The coupling strengths for π transitions are found by scaling (40) to be with currents of amplitude I (ω ↑↓ ±ω X/Y ) blue at frequencies ω ↑↓ ± (ω X/Y +δ X/Y ) that are required to drive Mølmer-Sørensen-type interactions [2]. For σ ± -transitions analogous relations hold involving the projections of µ d along (X ± iŶ )/ √ 2. We conclude that all interactions are similar in magnitude and that the current amplitudes can be used to tune the effective coupling strengths J X , J Y , and J Z of the Kitaev model simulator (36). To drive all bonds simultaneously, a total of five alternating currents at different frequencies are necessary; using the largest allowed value of 1 in (41), the maximum rms current each of the blue wires [and red wires, see (38)] has to sustain is I 2 blue (t) ≈ 5/2 × 0.36 A × [δ/(2π kHz)] 3/2 . We recall, however, that (40) and (42) depend very strongly on the vertical distance h w , and even a small reduction in h w can substantially reduce the required currents.
Expression (40) seems to suggest that for quantum simulators built with the principles described here, decreasing the physical size of the ion-trap lattice (as given by the length scale d) will strongly increase the simulation speed, which is given by the effective dynamics of the particular quantum simulator but ultimately proportional to the spin-spin coupling strengths. However, a careful analysis of assumptions and constraints, carried out below for interactions along Z but equally valid for all other effective interactions, disproves this observation. Firstly, if we assume that avoiding ion motion is a significant experimental constraint, a constant ratio J Z /( δ Z ) in (41) implies that the effective spin-spin coupling strength scales as for a given ion species and given electrode shapes. Secondly, assuming the currents to be limited by heat dissipation, an upper bound on I (ω Z ) blue must scale as d 3/2 . And lastly, lower bounds for the trap frequency can be found in two ways: (i) to meet our assumption of stiff trapping, i.e., ω 0Z ω Z , we requireω Z d −3/2 × |Q|/ √ 8π 0 M ; and (ii) to use the expansion (26) we require ω 0Z |δ Z |, which, combined with the scaling of (43) for δ Z and with the above current scaling, implies a scaling of d −5 for the lower bound of ω Z . Together, these bounds imply that the maximal achievable coupling strength (43) scales only as d −1/3 or even d 2 , depending on which of the frequency bounds is more stringent. Since current experimental setups are far from reaching these lowers bounds onω Z , miniaturization is expected for now to increase the coupling strength faster than these estimates; but the optimal dimension, where the ratio of simulation speed and heating rate (anomalous heating, scaling as d −4 [44,45]) is maximized, remains an open question.

Conclusions
We have discussed modifications to Coulomb potentials and interactions of trapped ions due to the presence of trap electrodes and cover planes. For plane geometries we have treated these modifications rigorously, using the method of image charges. We have found considerable deviations of the long-range behavior from that in free space when the relevant distances are of the order of ion-to-surface distances or larger. Moreover, we have developed a general approach to treating the effective spin-spin interactions of ions trapped in a multi-trap array in the stiff-trapping limit, where dipole-dipole interactions between nearest neighbors produce only small corrections to the bare normal modes of a given trap well. We have shown that effective coupling strengths, and therefore simulation timescales, are determined by the nearest-neighbor dipole-dipole couplings. As an illustration of the versatility and power of this stiff-trap-array approach, we have discussed a quantum simulation of the hexagonal Kitaev model. We have also addressed several practical challenges, including how the trap depth of the array may be improved so ions created from a thermal source with large kinetic energies can be trapped.