Quantum simulations with complex geometries and synthetic gauge fields in a trapped ion chain

In recent years, arrays of atomic ions in a linear RF trap have proven to be a particularly successful platform for quantum simulation. However, a wide range of quantum models and phenomena have, so far, remained beyond the reach of such simulators. In this work we introduce a technique that can substantially extend this reach using an external field gradient along the ion chain and a global, uniform driving field. The technique can be used to generate both static and time-varying synthetic gauge fields in a linear chain of trapped ions, and enables continuous simulation of a variety of coupling geometries and topologies, including periodic boundary conditions and high dimensional Hamiltonians. We describe the technique, derive the corresponding effective Hamiltonian, propose a number of variations, and discuss the possibility of scaling to quantum-advantage sized simulators. Additionally, we suggest several possible implementations and briefly examine two: the Aharonov-Bohm ring and the frustrated triangular ladder.


INTRODUCTION
Quantum simulators are highly controlled quantum machines with which it is possible to engineer and study complex quantum states and dynamics. Such machines, when large and accurate enough, are expected to elucidate the behaviour of quantum systems that defy analytical treatment and which are intractable for classical numerical simulations [1]. With the increasing sizes and abilities of quantum simulators and computers [2][3][4][5], quantum advantage in the context of quantum simulation may be within reach in the near future [6,7]. Of the diverse physical platforms used for quantum simulation, atomic ion chains in linear RF traps have proven particularly fertile by virtue of their long coherence times and high operation fidelity [8,9]. Using trapped ion quantum simulators, researchers have created and studied a wealth of quantum phenomena by applying both analog [3,[10][11][12][13][14][15][16][17][18][19][20][21] and digital [22][23][24][25][26] simulation techniques.
A principal feature of ion chain quantum simulators is the precisely controllable long-range coupling between ion-qubits, driven by an external field and mediated by the motional modes of the chain [27][28][29]. Representing a spin- 1 2 particle by two electronic energy levels in each ion, a uniform external driving field can induce an effective spin-spin interaction of the form [30]: where i, j denote the spin index, σ α i with α ∈ {x, y, z} are the standard Pauli operators acting on spin i, and J α ij represents the coupling matrix for the different Pauli axes. Often, quantum simulation experiments with trapped ions use a uniform bichromatic field to couple the spins through the transversal motional modes of the ion chain, generating a coupling matrix [9,13,30]: Despite the effectiveness of these tools, many territories remain uncharted for linear ion trap quantum simulators.
One outstanding challenge is that of simulating systems in more than a single spatial dimension. The successes of 1D ion trap quantum simulators calls for extending their scope to explore the richness of higher dimensional quantum systems. However, ion chains are open-ended and one dimensional, and the couplings that are induced by the simplest and most robust simulation techniques, expressed in Eq. (2), naturally reflect this geometry. In the past several years, new methods have been developed in order to enable more complex coupling geometries [31][32][33][34][35][36]. Despite these advancements, hardly any scalable simulations of high dimensional Hamiltonians have been shown in an ion chain.
Another tool that linear ion trap simulators currently lack, yet may aspire to, is the simulation of magnetic fluxes. Magnetic fluxes are a key ingredient in a range of quantum phenomena, with the iconic example being the quantum Hall effect [37,38]. Such fluxes can be expressed in the Hamiltonian by complex coupling terms, such as e iφ ψ † i ψ j , representing the presence of a gauge field potential which associates a phase to a directional propagation of an excitation along the lattice, also known as a Peierls phase. For this reason, the synthetic or artificial generation of gauge field interaction terms has emerged in recent years as one of the most prolific tools of neutral atom quantum simulation [39][40][41][42][43][44][45][46][47], and might similarly present new opportunities in trapped ion quantum simulators.
In principle, making use of the universal gate set already available in ion trap quantum computers [23,48,49], one can digitally simulate any quantum system by breaking down the dynamics into a series of simpler operations [8,23,50]; such a simulation can include all features discussed above. However, in practice, the engineering cost of a universal gate set is high and the decomposition of target models may be unwieldy and can incur a high fidelity cost. While in the long term digital simulations may benefit from fault-tolerant quantum error correction, the necessary qubit array sizes and operation fidelities to reach this threshold are far beyond current capabilities. Hence, analog quantum simulations, in which the target Hamiltonian is continuously implemented, arguably offer more promising prospects for near and mid-term quantum simulation [51]. This motivates expanding the range of models that are directly simulatable with trapped ions.
In this manuscript we introduce a scalable and experimentally simple technique that can be used to simulate a large range of spin Hamiltonians on an ion chain. This technique improves on the standard schemes in two significant ways: through generation of complex coupling geometries, including high dimensional Hamiltonians and closed boundary conditions; and by an introduction of both a static and time-dependent Peierls phase, effectively generating a synthetic gauge field. Crucially, the technique can be performed with a uniform intensity global driving beam, with no dynamical control, and with fields that are independent of the number of ions. The technique requires an addition of an external gradient field along the chain and the use of a multitone driving field. Gradient fields have been most prominently used in ion chains in quantum processing architectures where a strong spectral separation enables both individual addressing of the ions (as in NMR) as well as driving entangling gates using long-wavelength fields [52][53][54]. The technique outlined in this manuscript can be understood as an extension or variation of the NMR-inspired scheme, spectrally resolving coupling terms rather than individual subsystems. A similar proposal for employing gradient fields to simulate higher dimensional systems was recently put forth by Rajabi et al. [32], requiring the additional use of dynamical techniques.
We note exciting proposals for simulation of static [55,56] and dynamic [34,57] gauge fields as well as high dimensional Hamiltonians [31,34] in a trapped ion chain using either dynamical techniques [55,56], individual addressing of all ions [31,34], or additional energy levels [57]. Furthermore, proposals have been put forth for implementing synthetic gauge fields for the motional, rather than electronic, degrees of freedom of the ion chain [58][59][60].
The manuscript is ordered as follows: we first present our main result; we then describe the technique's principle of operation and derive the coupling Hamiltonian; we use the Aharonov-Bohm ring in order Here we only plot a pair of blue sidebands for illustration purposes, where the corresponding red sideband pair is not shown. The controlled resonances can be used to tailor different coupling geometries, such as 1D rings (c). The phases of the driving pairs can be chosen to generate synthetic gauge fields, representing magnetic fluxes threading the lattice. In this example, as a result, an excitation on the lattice will be driven around the ring, in a direction and velocity dictated by the flux (d). Here the flux is Φ = 3π/4. The probability Pe for a local spin excitation is color coded.
to exemplify the salient features of our technique; we suggest a number of simulatable Hamiltonians of interest; and finally, we discuss the challenges posed for realizing the technique on large quantum simulators. In the appendix we describe variations that can ease implementation and further increase the range of target models.

MAIN RESULTS
Our main result is a simple recipe for generating a class of Hamiltonians using trapped ions. The class is described by the following formula: where σ + i (σ − i ) denotes the raising (lowering) Pauli operator on ion i, and Ω n , φ n and δ n are tunable parameters corresponding to the coupling strengths, static phases and time-dependent phases, respectively, of an n-neighbor hopping interaction. As we show below, the Hamiltonian in Eq. (3) can be used to implement spin Hamiltonians on various geometries, and can furthermore manifest static and time-dependent gauge fields.
The method relies on the application of a static external field gradient (e.g. a spatially varying magnetic or light shift field) for shifting the atomic energy levels along the ion chain, applied together with a corresponding global uniform driving field. The field gradient collapses the translational symmetry of the ion chain, effectively suppressing the standard coupling form of Eq. (2). However, because the gradient is spatially uniform, the driving field can be used to selectively reinstate the translational symmetry of the spin-spin interaction. This is done in a controlled manner by bridging the resonance difference between equally-separated ion-pairs using the frequency difference between pairs of bichromatic fields.
Furthermore, the breaking of spatial symmetry differentiates the interaction of an ion with its neighbors to the left and right, which gives rise to a gauge-field-like phase and an effective breaking of time-reversal symmetry.
The tools needed to implement (3) are standard in trapped ion experiments: for every nonzero Ω n , representing an n-neighbor interaction, a four-tone field is added. The corresponding Ω n , φ n and δ n are set by the field's amplitudes, phases and frequencies. The driving field is activated using a single uniform-intensity beam. The magnetic field gradient can be modest, on the order of 10 G/cm. Additionally, the interaction in Eq. (3) is excitation-number-preserving implying robustness to global dephasing noise.
A wide spectrum of quantum phenomena can be accessed using this method. For example, by choosing Ω 1 = Ω N −1 = Ω and φ 1 = −φ N −1 = 2πΦ/N (and nulling all other parameters) we arrive at a lattice ring Hamiltonian with a hopping term: where boundary conditions are periodic. The phase Φ corresponds to the Aharonov-Bohm phase acquired by an electric charge encircling a ring penetrated by a magnetic flux (note that a 1D spin system may always be described as a fermionic system, using the Jordan-Wigner transformation) [61]. Accordingly, an excitation will travel clockwise or counter-clockwise on the ring, generating a persistent current [62]. While this model is easy to solve, it clearly showcases the main tools of the proposed technique. We analyze the model in more detail below.
Figure (1) highlights our method's main principles of operation. Taking the 1D ring in a N = 5 ion chain as an example, a pair of driving fields (a) bridge the energy difference between the ion created by the external gradient, i.e ∆ for neighboring ions and 4∆ for the edge ions (b). These driving fields form tailored couplings between the ions, and are here used to generate a 5-site ring penetrated by a magnetic flux Φ (c). Accordingly, a simulation of the ion chain's evolution shows an excitation travelling around the ring (d).
Using these principles, Hamiltonians of the form expressed in Eq. (3) can be generated. Significantly, a wide variety of coupling geometries are reachable. We illustrate some possible coupling geometries in Figure  2. Besides a ring (a), these include triangular ladders (b); 2d rectangular lattices (c) (which can be closed onto a cylinder, not shown); a Möbius-strip ladder (d); a helical lattice on a cylinder (e); and a torus (f). These lattices can be threaded by a variety of magnetic fluxes, as illustrated for the torus (f). While the 1D ring can be mapped to a system of free noninteracting fermions and is thus simply solvable, all other models shown here are expected to show complex behavior and can be difficult to solve.

PHYSICAL PICTURE
Ions in a Paul trap are frequently modeled as two level spins with a set of harmonic modes, where the former corresponds to the ions' electronic degrees of freedom, and the latter to the motional normal-modes of the ion chain. External electromagnetic fields can couple to spin and motional degrees of freedom and, with proper tuning, can be used to engineer effective interactions between the spins of different ions via mediation by the motional modes.
In the Mølmer-Sørensen (MS) interaction [27,28], the external field is bichromatic and tuned to frequencies ω ± = ω 0 ± (ν + ξ), with ω 0 the single qubit energy separation, ν the frequency of a normal mode of motion of the ion-chain, and ξ a constant detuning which together with the field intensity determines the interaction rate. The tone ω + (ω − ) mediates interactions via the blue (red) motional sideband, i.e it employs transitions which excite the ion's electronic degree of freedom while adding (removing) a phonon of the motional normal-mode. While only two driving tones are used, this interaction couples any two ions in the chain in four different "pathways", as is shown in Fig. 1 of Ref. [28].
When coupled to the center-of-mass (COM) mode of the ion chain, this driving field induces an effective σ x σ x interaction through a two photon process, which equally couples all of the ions-pairs in the ion chain. This interaction can be decomposed to two contributions: a pair creation/annihilation term, σ + σ + +h.c, which drives a b d c e f Figure 2. Implementation of various geometries using appropriate hopping interactions.
(a) nearest-neighbour and N − 1 neighbor interactions generate a ring; (b) nearest-neighbour and next-nearest-neighbour interactions generate a triangular ladder; (c) nearest-neighbour and W -neighbor interactions generate a rectangular lattice when used with spacer ions (see discussion below), and with the addition of N/W interactions generate a rectangular lattice on a cylinder; (d) nearest-neighbour, N/2 and N − 1 interactions generate a Möbius ladder; (e) nearest-neighbour and W interactions generate a helical lattice on a cylinder; (f) by adding N − 1 and N/W interactions to the helical lattice, the cylinder is closed onto a torus. Controlling the phases of these interactions results in synthetic gauge fields representing fluxes threading these geometries. For instance, in the torus (f) an external axial flux (green), or within the torus, as in an anapole moment (red), can be produced.
Here we are interested in eliminating the pair creation/annihilation term while retaining the hopping term, and furthermore shaping its coupling matrix. The first goal is achieved by detuning the driving field frequencies from resonance with the two-photon transition, i.e by modifying the bichromatic drive frequencies to ω ± = ω 0 + ± (ν + ξ), shifting the pair creation/annihilation term 2 off-resonance. The second goal requires a more elaborate approach. As the σ + σ − + h.c term is mediated by an excitation/de-excitation pair of identical photons, it resonantly couples only states that are degenerate under H 0 ; this implies a coupling between all equal-excitation states. However, the hopping term can also be controllably suppressed by lifting the equal-excitation degeneracy [63]. A controlled suppression of the interaction will then allow for selectively reinstating resonant conditions through a modulation of the driving field.
To do so, an external (e.g. magnetic) field gradient is added along the ion-chain, such that the transition frequency between adjacent ions differs by ∆. In order to selectively couple ions which are n sites apart we drive the ions with four frequencies, composed of the frequency pairs ω b,± = ω 0 + (ν + ξ b ) ± ∆n 2 and ω r,± = ω 0 − (ν + ξ r ) ± ∆n 2 . The pair ω b,± couples an n-site hop resonantly, mediated by the blue sideband. Similarly, the pair ω r,± couples the same hop, mediated by the red sideband. As in the MS interaction, both sidebands are employed in order to mitigate temperature-dependent effects. In order to keep the pair creation/annihilation term non-resonant we use ξ b = ξ + and ξ r = ξ − . Figure 3 illustrates the transition from all-to-all coupling in the absence of a gradient field (a) to a complete suppression of coupling due to the gradient (b) and the selective resurrection of coupling by introducing the resonant sideband modulation (c,d). DERIVATION We outline the derivation of Eq. (3). We focus only on a single H n term, and later comment on the generalization to a summation of these terms. In the absence of an external driving field, the Hamiltonian of N trapped ions in a magnetic field gradient is given by: Here ∆ is the transition energy difference between adjacent ions due to the field gradient, and ν l is the frequency of the l-th normal-mode with the annihilation operator a l . In our derivations below we will assume coupling to a single motional mode, the COM mode (with a frequency ν = ν 1 and a Lamb-Dicke parameter η = η 1 ), which couples equally to all ions in the ion-chain. This assumption can be relaxed [33], as will be shown in later discussions. For each n we apply a four-tone field composed of the frequency pairs ω b,± = ω 0 + (ν + ξ b ) ± ∆n 2 and ω r,± = ω 0 − (ν + ξ r ) ± ∆n 2 , Rabi frequencies Ω b and Ω r , and phases φ b,± = ±φ/2 and φ r,± = ± (φ + π) /2. In addition ξ b and ξ r are chosen such that pair creation/annihilation transitions between all ion-pairs are detuned from resonance.
Specifically, = (ξ b − ξ r ) /2 is chosen such that the detuning from any pair creation/annihilation resonances, i.e the transition frequency of any two states which differ by two excitations, is large compared to the effective coupling η 2 Ω r Ω b /ξ [28,63]. Thus, the red and blue sideband pairs contribute to the evolution independently.
We first focus on the interaction mediated by the blue sideband, which is due to the pair ω b,± . In a frame rotating with respect to H 0 above, this interaction is described by, This expression is valid in terms of a rotating wave approximation by assuming that Ω b /ν and Ω b /ω 0 are small, and in leading order in η.
The Hamiltonian in Eq. (6) cannot be solved analytically. However, its resulting evolution operator, U (t), can be approximated by using the leading terms in a Magnus expansion [64,65], U = exp − i χ (t) = exp − i n χ n (t) . Since χ n ∝ ηΩ b ξ n , we are satisfied with terminating the expansion at the second order, which (as we show below) provides the leading order resonant terms.
The evolution due to χ is stroboscopic with a fundamental period T , i.e that χ (kT ) = kT H eff , with k ∈ Z and H eff an effective time-independent Hamiltonian. In the limit T → 0 the derived Hamiltonian approaches the target Hamiltonian at all times.
In first order the expansion reads, By choosing T ∆ = 4πm and T ξ b = 2πM b , with m, M b ∈ Z, we arrive at χ 1 (T ) = 0, trivially satisfying the stroboscopic condition.
In the next order we decompose χ 2 to a hopping term and a rotation around the z-axis. The expansion is then given as, which can be translated to an effective Hamiltonian, with corrections that scale as N ∆ ξ b 2 . The first two terms of this Hamiltonian are homogeneous. The first term represents the desired hopping interaction (hence the subscript h), while the second is effectively equivalent to a temperature-dependent global magnetic field in the z-direction (hence the subscript z), with Ω n,b = Since H h is excitation preserving, by initializing the system to an eigenstate of k σ z k , i.e. to a state with a well defined number of excitations, H z is reduced to a global phase and can be ignored. For these initial states a two-tone driving field suffices.
Both H h and H z are proportional to Ω n,b , however the former also depends on the phase φ. This enables the use of the red sideband pair, ω r,± , in order to eliminate H z entirely while maintaining H h . To this end, we choose ξ r and Ω r such that Ω n,r = −Ω n,b , and φ r,± = ± φ+π 2 . Thus the combination of the two pairs yields H h → 2H h and H z → 0.
The two latter terms in Eq. (9), corresponding to the two terms in Eq. (10), are These terms are gradient inhomogeneous terms, i.e they are not tranlationally invariant, and vanish in the limit ∆ → 0. By setting ξ b N ∆, H ∇ h and H ∇ z become negligible and we obtain a homogeneous effective Hamiltonian. Furthermore, by using the red sideband pair as described above H ∇ h is eliminated entirely as well. The remaining H ∇ z is more difficult to eliminate in the non-adiabatic regime. If the ion chain is cooled to the ground state, its contribution to the effective Hamiltonian is simplified to H k , which has the same form as that of the external gradient field in (5). Hence, by making a small correction to the addressing field frequencies the ground-state contribution of H ∇ z is eliminated. Thus, at the appropriate limits, the four-tone frequency drive yields the effective Hamiltonian, The interaction may be made time-dependent by detuning the two-photon transition, ω b,± → ω b,± ± δ/2, with δ ∆ (similarly for ω r,± ). We obtain an off-resonant coupling, which manifests in a time variation of the hopping phase: σ + k+n σ − k e iφ → σ + k+n σ − k e i(φ+δt) . With this transformation we obtain the general form of H n from Eq. (3).
We further note that more generally the Peierls phase can be changed in time in whatever way one wishes by changing the appropriate driving field phases, as long as all spectral component of the phase dynamics are much smaller than ∆. Constant detuning, implying a linear change in time of the phase, is a specific instance of such phase dynamics. As an example, by periodically modulating the phase, the presence of AC magnetic fluxes can be realized.

SOME TARGET MODELS
In this section we briefly explore two models that can be simulated with our method and suggest several additional models. In order to clearly display the salient features of the technique, we first focus on the analytically solvable 1D Aharonov-Bohm ring. We then discuss a triangular spin ladder model exhibiting geometric frustration. Despite the simplicity of this model, it gives rise to a rich phase diagram and interesting physical phenomena.

The discrete 1D Aharonov-Bohm ring
For N ions, turning on the n = 1 and n = N − 1 interaction terms in Eq. (3) generates a 1D ring lattice. Adding phase terms φ 1 = −φ N −1 = 2πΦ/N simulates a magnetic flux Φ penetrating the ring, as is shown in Fig. 1, with the Hamiltonian H r expressed in Eq. (4). As the Hamiltonian is excitation preserving, inside an excitation eigenspace H r can be mapped directly via the Jordan-Wigner transformation [66] to fermions on a ring lattice threaded by a magnetic field: with periodic (antiperiodic) boundary conditions for an odd (even) number of excitations. H f is the spinless fermion 1D tight-binding model [67]. The magnetic flux threading the ring may give rise to a persistent current due to the Aharonov-Bohm effect, which survives even in the presence of impurities in the chain [68,69]. In order to observe this effect, one can prepare an initial state with a position occupation distribution that will rotate around the ring without diffusing. The singly-excited subspace of H f is spanned by plane-wave eigenstates |k = 1 √ N n e ink/N ψ † n |0 with energies E k = 2 cos 2π N (k + Φ) . While each of these waves uniformly occupies all ions along the chain, we can initialize the system in a wave packet state with a position-dependent occupation: |ψ W.P (t = 0) = |k + e iϕ |k − 1 / √ 2 [70]. In this state the probability to occupy the n-th site is Equation (14) shows that |ψ W.P is a wave-packet, with ϕ determining its position on the ring. The wave-packet's evolution can be described by the evolution of ϕ. The state will evolve according to That is, the packet rotates around the ring at a constant, flux-dependent, angular velocity. As expected, at the large N limit, we obtain v (Φ) ∝ Φ + k. The Aharonov Bohm ring can be used to observe Bloch oscillations. For particles in a 1D periodic structure, the addition of a constant uniform force generates an oscillatory motion rather than unidirectional acceleration [71]. In a 1D ring such a force can be created by threading the ring with a time-dependent magnetic flux [68]. An excitation, rather than encircling the ring with a constant acceleration, will oscillate locally. The effect can be naturally incorporated using our technique, taking advantage of the ability to generate a time-varying synthetic gauge field using off-resonant driving pairs within the coupling scheme of the AB ring, as previously described.

Triangular Spin Ladder
The Aharonov-Bohm ring can be mapped onto a free fermion model via the Jordan-Wigner transformation, and is thus easily solvable. However, spin interactions beyond nearest-neighbor can only be mapped onto interacting fermion models, and accordingly generate complex dynamics and phases which are often challenging for classical computation techniques. As an example, we briefly discuss the triangular ladder, a simple Hamiltonian that can be easily implemented using our technique, but which nonetheless manifests complex behavior.
Activation of the n = 1 and n = 2 terms in Eq. (3) generates a nearest neighbor (nn) and next-nearest neighbor (nnn) interaction spin Hamiltonian: where boundary conditions are open. Such Hamiltonians can be graphically represented by triangular ladders in which rungs and rails represent nn and nnn interactions accordingly, as pictured in Figure 5. Due to the competition between nn and nnn terms, geometrically viewed as the competition of interactions inside each triangle, models of this sort are frustrated and thus give rise to a relatively rich phase diagram [72][73][74]. H tl is gauge invariant under the transformation φ 1 → φ 1 + ϕ, φ 2 → φ 2 + 2ϕ for any ϕ; this is equal to the gauge transformation σ + k → e ikϕ σ + k . For trivial interaction phases φ 2 = 2φ 1 and antiferromagnetic nnn interactions (J 2 > 0), H tl represents the one-dimensional frustrated XY chain model for spin-1 2 . This model is a paradigmatic example of frustration [73]. It supports a variety of phases, notably including an exotic chiral-ordered phase [75][76][77][78]. The system phase depends on j = J 2 /J 1 , which can be fully controlled using the techniques outlined in this manuscript. At j = 1 2 , also known as the Majumdar-Ghosh point, the ground state is an exactly solvable dimerized state [74]. The control and flexibility of trapped ion systems may enable generation of these unique phases and direct measurement of their order parameters [79], their entanglement properties [80] and their excitation dynamics [16].
By choosing nontrivial values for φ 1 , φ 2 an additional synthetic gauge field representing a staggered flux is added to the Hamiltonian. Each triangle is pierced by a gauge-invariant magnetic flux ±Φ = φ 2 − 2φ 1 , with the flux alternating signs between neighbouring plaquettes. Similar triangular ladder models with synthetic flux fields have been suggested and implemented in neutral atom systems [81][82][83][84]. Figure 5 illustrates the connectivity and staggered flux for this model. Figure 5. The triangular spin ladder. By turning on both nn and nnn couplings, an effective triangular or zigzag ladder coupling geometry is generated. The triangular layout can easily lead to frustration, which can be understood as a competition of the order imposed by the different coupling terms. Corresponding coupling phases φ1 and φ2 create a synthetic gauge field representing a magnetic flux that alternates in sign between plaquettes, with a gauge-invariant phase Φ = φ2 − 2φ1.

Rectangular lattice
In the 2D examples discussed in this manuscript so far, lattices were either triangular or helical, and not rectangular. This is due to the fact that a strictly 2D rectangular lattice cannot be reduced to the form given by Eq. (3); placing qubits on the lattice, the nearest-neighbor coupling scheme would create a link between the last qubit in row k and the first qubit in row k + 1, violating the lattice geometry. This can be remedied by interrupting the nearest-neighbor interaction chains, represented by the rows of the lattice, through insertion of an auxiliary passive ion. The auxiliary ion acts as a spacer, generating an effective jump in the gradient for the active ions. This ion can be of a different isotope or species, but is most easily chosen to be an identical ion that is either strongly light shifted by an individual addressing beam or prepared in any state outside the qubit subspace.
In the simplest example, a single spacer ion in the middle of a chain of 2N + 1 ions, along with interaction terms H 1 and H N +1 , would generate the rectangular spin ladder, shown in Figure 6. With the addition of more spacer ions, more rows could be added to this array, effectively creating a complete rectangular lattice.
This lattice can then be curled into a cylinder with the activation of an additional term. For N = h · (w + 1) ions, where h, w ∈ Z, activating the nn, w + 1 and N − 1 − (w + 1) terms generates an h × w rectangular lattice a b Figure 6. Spacers and rectangular ladders. (a) Placing the spins on a rectangular chain and activating the n = 1 and n = m terms generates a rectangular lattice with additional unwanted cross-row terms. (b) By adding spacer ions, initialized in a state which is uncoupled to the spin dynamics, and by replacing the n = m term with a n = m + 1 term, the unwanted coupling is corrected, giving an exact rectangular coupling geometry.
on a cylinder. The curled dimension can be threaded by a flux, determined by phases of the non-nn terms.

Additional geometries
There are a number of other geometries which can be generated in a straightforward manner using our technique. We briefly mention several more examples: the Möbius ring, the cylindrical helix, and the torus, all illustrated in Fig. 2.
Turning on the terms H 1 and H N/2 generates a 2×N/2 rectangular lattice with an additional connection of site n−1 to site n, as is shown in Fig. 6a for n = 3. Activating in addition the term H N −1 forms a Möbius ring, shown in Fig. 2d. Such a system may be used to study topological effects in non-trivial geometries [85][86][87][88].
For N = w · h where w, h ∈ Z, turning on the terms H 1 and H w in Eq. (3) generates the cylindrical helix in Fig. 2e, with w sites per loop and height h. Furthermore, adding terms Ω N −1 and Ω N −w induces periodic boundary conditions, resulting in the torus seen in Fig 2f. Setting the phases φ 1 = −φ N −1 and φ w = −φ N −w gives rise to two independent fluxes penetrating the torus, Φ 1 = Lφ w (green arrow) and Φ 2 = W φ 1 − φ w (red arrow). Such a system may be used for the study of the quantum Hall effect in the thin torus limit [89,90].
In Appendix A we show a simple resource efficient implementation of our method, that can be used for the realization of some of the models above.
The ideas presented here may be taken even further by adopting the powerful neutral atoms quantum simulation concept of synthetic dimensions [47,91,92]. In neutral atom quantum simulators, extraneous internal degrees of freedom of the atom are used to represent additional lattice sites. For instance, a 1D system of atoms can be used to represent a 2D lattice, where the supplementary dimension is embodied by additional internal states of each atom. In such a case, engineered spin-orbit coupling can be used to drive a synthetic gauge field. In trapped ions, the additional Zeeman or hyperfine states may present a similar possibility, further extending ion chain simulations to an additional dimension.

Spatially varying potentials
In ion chains, it is possible to generate site-dependent energy shifts, H v = 1 2 k V k σ z k , by using individual addressing beams, nonuniform global beams, magnetic fields or other spatially varying fields.
Under the condition V k ∆, H v can be applied in parallel to our technique. As H v commutes with H 0 , it can be simply added to the interaction picture Hamiltonian, V I → V I + H v . In this instance, H v can be interpreted as a spatially varying potential on the lattice. By choosing random site-dependent shifts, disorder is added to the system. Disorder can give rise to localization [93], which can thus be studied in the variety of contexts presented in the paper.

SCALING UP
While implementation of the technique we propose in small scale simulators should be straightforward, approaching larger, quantum-advantage (NISQ) sized simulators [51] seems feasible yet more demanding. Here we discuss possible challenges in the implementation of the technique in large quantum simulators.
The technique, as presented, calls for constant differences in resonance frequencies between neighboring ions. This can be exactly achieved with a spatially uniform gradient only to the extent that the ions are spaced equidistantly. However, many linear ion traps currently use a harmonic axial potential, in which ions are not equidistant [94]. Higher order components of the external field can be engineered to meet this issue, albeit with increasing experimental complexity. Nevertheless, several groups have constructed -or are in the process of constructing -anharmonic traps, with the stated purpose of trapping ions equidistantly [34,[95][96][97]. Anharmonic traps are likely to become more useful in future ion trap simulators and computers, due to their advantage in maintaining high inter-ion spacing, critical for preventing cross-talk in addressing and detection; their resistance to transitions from linear to zig-zag crystal configurations; and their suppression of inhomogeneous quadrupole shifts. Our proposal is best suited for such traps.
Another hurdle for large scale implementation could be the inverse relation between coupling strength and number of ions, keeping driving field intensity constant. We define adiabaticity parameters α = ξ/ (∆N ) and β = ∆/ η 2 Ω 2 0 2ξ . In the low α limit the inhomogeneous contribution of the gradient field to the Hamiltonian becomes prominent (although strongly suppressed by the double sideband frequency configuration), and for α < 1/2 some ions in the chain may even be resonantly driven on the sideband transition. Similarly, for small values of β, non-resonant Hamiltonian terms, which are ideally completely suppressed, can play a significant role in the dynamics.
Keeping these adiabaticity parameters constant, the effective coupling strength is 2αβN . Since we assume coupling only to the COM mode, the Lamb-Dicke parameter also decreases with ion number: η = η1 √ N , where η 1 is the single-ion Lamb-Dicke parameter. Hence, overall, the coupling strength decreases as 1 N , as opposed to the normal COM MS degradation of 1 √ N . In essence, the additional penalty emerges from the requirement to preserve spectral spacing in presence of the gradient field. The degradation can be remedied by increasing the field intensity or by coupling to many modes rather than just the COM, which is not only convenient but necessary when using radial modes in large simulators.
The radial modes of linear ion traps bunch up when adding more ions, and consequently the radial COM mode cannot be spectrally resolved. Using these modes as interaction mediators will thus inevitably generate coupling to a multitude of modes. Even so, the use of radial modes can be advantageous and is compatible with our proposal, up to a modification of the effective coupling, in analogy with [30] (see Appendix B for more details).
In contrast, the axial COM mode remains spectrally separated from all other modes independently of the number of ions, and thus is in this sense ideal for the proposed technique. Nevertheless, axial modes come with their own disadvantages. Working with large ion chain requires very low axial trap frequencies in order to avoid the crystal zig-zag transition, or buckling, in the middle of the chain. Low axial frequencies lead to higher carrier coupling, temperatures, and heating rates, limiting simulation fidelity. These problems can be strongly mitigated by using anharmonic trapping potentials, which are beneficial for our proposal as stated above.

SUMMARY
In conclusion, we have introduced and explored a technique that realizes a variety of spin Hamiltonians in trapped ion chains. The range of implementable Hamiltonians includes spin lattices with dimension larger than one, closed boundary conditions, rectangular and triangular lattices, and full control of nearest and next-nearest neighbor couplings. Furthermore, the technique provides a means to realize static and time-varying synthetic gauge fields in ion chains. This is done using a global driving field of uniform intensity, with no need for individual addressing or for dynamic control, and a static external field gradient along the ion chain. The advanced tools that have already been developed for trapped ion chains, including preparation of highly entangled states [98] or measurement of observables of interest such as entanglement entropy [80], encourages us to believe that ion chain analog simulation can break new ground in territories which have been considered outside its purview, such as simulation of 2D topological phenomena. The tools we present here are a step towards this direction. generalization of Eq. (8), it is given by, χ 2,j = χ 2,j;h + χ 2,j;z χ 2,j;h = iT 2 Ω 2 b k B k,k+n σ + k+n σ − k e iφ + h.c + O α −2 χ 2,j;z = 2iT 2 Ω 2 b j,k B k,k+n σ z k a † j a j + where we defined B i,k = j η j,i η j,k /2ξ b,j . That is, we obtain the same effective Hamiltonians as in Eq. (10) but with the normalized spin-spin coupling B i,k , which comes about due to contributions from all of the normal-modes, in analogy to [30].
Since the radial modes are bunched, typically |ξ b,j − ξ b,k | |ν j − ν k |. Consequently the red sideband can be used in order to eliminate H z , as described in the single normal-mode case above.
The term χ 2,j,j couples between mode j and mode j . It only contains operators of the form σ z k a † j a j , and its conjugate; therefore, this term generate no spin-hopping. By using the red sideband frequency pair, the leading order contribution of this term scales as α −2 and is therefore neligible in the large α limit.
We thus conclude that our method, in its continuous instance, is compatible with coupling through the radial motional modes, with the appropriate modification of coupling strengths.