Hexagonal Plaquette Spin-spin Interactions and Quantum Magnetism in a Two-dimensional Ion Crystal

We propose a trapped ion scheme en route to realize spin Hamiltonians on a Kagome lattice which, at low energies, are described by emergent Z2 gauge fields, and support a topological quantum spin liquid ground state. The enabling element in our scheme is the hexagonal plaquette spin-spin interactions in a 2D ion crystal. For this, the phonon-mode spectrum of the crystal is engineered by standing-wave optical potentials or by using Rydberg excited ions, thus generating localized phonon-modes around a hexagon of ions selected out of the entire two-dimensional crystal. These tailored modes can mediate spin-spin interactions between ion-qubits on a hexagonal plaquette when subject to state-dependent optical dipole forces. We discuss how these interactions can be employed to emulate a generalized Balents-Fisher-Girvin model in minimal instances of one and two plaquettes. This model is an archetypical Hamiltonian in which gauge fields are the emergent degrees of freedom on top of the classical ground state manifold. Under realistic situations, we show the emergence of a discrete Gauss's law as well as the dynamics of a deconfined charge excitation on a gauge-invariant background using the two-plaquettes trapped ions spin-system. The proposed scheme in principle allows further scaling in a future trapped ion quantum simulator, and we conclude that our work will pave the way towards the simulation of emergent gauge theories and quantum spin liquids in trapped ion systems.


Introduction
Topological quantum spin liquids are fascinating states of matter supporting topological order and exotic excitations with fractional statistics [1]. Recently, there has been intense theoretical activity aimed at underpinning their properties, and identifying microscopic Hamiltonians which could support such states. This activity is mostly focused on the weak and strong insulator scenarios [1]: in the first case, the focus is on SU(2) invariant spin Hamiltonian on frustrated geometries [2], which could potentially be relevant for a series of solid state materials. In the strong insulator scenario, one is instead interested in Hamiltonians displaying exotic constrained dynamics [1, 3,4]. This second route has been particularly fruitful at the theory level, as the corresponding model Hamiltonians are amenable to powerful analytical techniques and, most importantly, provide a very clean interpretation of quantum spin liquids in a gauge theory framework [1,5]. However, due to the exotic form of the basic constraints necessary to realize a strong insulator dynamics, many paradigmatic Hamiltonians in this class are still thought of as idealized dynamics without any direct physical counterpart.
Here, we show how a minimal instance of an archetypical example of spin Hamiltonian supporting a topological quantum spin liquid ground state, the Balents-Fisher-Girvin (BFG) model [6], can be realized in 2D laboratory trapped ion systems.
The key element is to engineer constrained dynamics in a controlled fashion, and in particular, to realize exotic 'plaquette' constraints typical of spin Hamiltonians in the block for quantum simulators aimed at studying the fundamental nature of frustrated quantum magnetism with emergent gauge fields [2,4,37] and may open up new directions in ion-based quantum simulators. By pinning multiple ions, more individual plaquette interactions can be engineered thus allowing for scaling up the quantum simulator. In the present paper, we discuss the minimal instance of a two-plaquette implementation, which should lie within current experimental reach. Note that our current studies also fit well into the recent trends of cold atom-investigations on synthetic gauge fields [38], in particular on dynamical gauge fields [39,40,41,42,43,44,45].
The ion pinning can be accomplished by a standing wave light field focused onto the ion to be pinned. Such optical potentials can reach curvatures corresponding to trapfrequencies in the MHz range. The proposed setup resembles those encountered in recent studies involving ions interacting with strong light-fields, such as in anomalous diffusion of an ion [46], preparing nonclassical motional states [47], and, in particular, modifying the trapping potential locally in a rf-trapped ion crystal for quantum simulations [48,49]. Additionally, optical trapping of an ion has been shown experimentally using either a single beam dipole trap [50,51,52] or an optical lattice [53,54,55,56]. An interesting alternative option would be to dress the ion of interest with a Rydberg state, such that the dipole moment induced by the ionic trapping field causes a change in local trapping frequency for the ion [57,58]. Here, no spatial variation of the pinning laser would be required and additional microwave fields could be employed for fine-tuning [59,60].
The paper is organized as follows: in section 2 we discuss the physical setup of the ion-laser system and the governing atom-light interaction Hamiltonian. The spin models are derived from the interaction Hamiltonians, we show how to generate spinspin hexagonal interactions in a single plaquette using N = 7 and two plaquettes using N = 19 ion crystals. In section 3 we discuss the ground state properties and magnetization dynamics for the single and double plaquettes spin-systems, and the results are compared with that of the original BFG model. Realistic numbers for the laser parameters for a calcium ion-setup are given in section 4. Finally, we summarize the paper with an outlook in section 5.

Spin-spin interactions in a two-dimensional ion crystal
A sketch of our proposed experiment is shown in figure 1. The system consists of a selfassembled 2D ion crystal confined in the xy plane by what we, for simplicity, assume to be a radially symmetric trap with trap-frequencies: ω x,y ω z . The internal level structure of the ions consists of three long-lived low-lying atomic states, these can be either three Zeeman states in the ground-state manifold or a combination of metastable and ground states. The states |2 ≡|↓ and |3 ≡|↑ encode the spin-half states of the BFG model. The ion crystal supports quantized phonon modes that are used to transmit The schematic picture of the ion-laser system for generating hexagonal-plaquette interactions. A pair of counter-propagating Raman fields with wave vectors (frequencies) k 1ẑ (ω 1 ) and k 2ẑ (ω 2 ) are used for inducing spin-spin couplings. OL1 and OL2 (yellow in colour) are two optical lattices focused on ions 1 and 2 to modify their transversal (along the z axis) trapping frequencies locally. (b) The level scheme: the internal level structure consists of three long-lived low-lying states and a manifold of excited states. The Raman lasers (blue and red arrows), detuned by δ m with respect to the frequency of a motional mode, generate state dependent optical shifts in the spin states |2 ≡ | ↓ and |3 ≡ | ↑ with an energy separation of ω ↓↑ , which simulate the spin-spin interactions of the type S i z ⊗ S j z between different ions i and j. The ions 1 and 2 prepared in state |1 experience an additional optical lattice that is created by very far detuned laser fields (yellow arrow) from the excited state |e with a detuning ∆ OL . The lattices modify their local trapping frequencies along the z axis. Additional laser fields, used for generating quantum dynamics in the spin system, are not shown in the figure for the sake of clarity.
spin-spin interactions between the ions via a Raman pair of counter-propagating laser fields as shown in figure 1b. Ions sitting in the middle of an hexagonal sublattice (see, e.g., ions 1 and 2 in Fig. 1) are prepared in the state |1 and strongly pinned by a standing wave laser field that is far detuned from any excited states {|e } to locally modify the phonon spectrum.
The details of laser mediated spin-spin interactions in trapped ion systems have been extensively described elsewhere [11] and are briefly discussed in Appendix A. Two counter propagating laser fields are far detuned from the transitions |↓ ↔ |e and |↑ ↔ |e . When the Rabi frequencies of these lasers are much smaller than the detunings, the excited state |e can be adiabatically eliminated. To engineer spinspin interactions of the form S z ⊗ S z , the frequency difference of the laser beams ω I is tuned close to the phonon frequencies ω m , where m denotes the particular phonon mode and δ m = ω I − ω m . Throughout the paper we assume to be in the Lamb-Dicke regime, η i m n m + 1 1, where η i m = q m0 b i m k I are the Lamb-Dicke parameters, n m the number of phonons in mode m, q m0 = /2M ω m with M the mass of an ion and b m = {b i m } are the (normalised) phonon eigenvectors. As long as η m Ω I (n) δ m , with the two photon Rabi frequency Ω I (n) = Ω 1 (n)Ω 2 (n)/∆ for each laser and state n =↑, ↓ and ∆ 1 ≈ ∆ 2 Ω 1,2 (n) the detuning for each laser, the Raman laser only transiently excites phonons and we can integrate out the phonon degrees of freedom altogether to obtain effective spin-spin interactions between ion i and j of the form with Ω i I = Ω i I (↓) − Ω i I (↑). The strength and range of J ij z are determined by the Rabi frequencies Ω i I , the detuning δ m and the amplitude of oscillation of each ion in the m th mode through the Lamb-Dicke parameters η i m . Note that we only assumed to be in the Lamb-Dicke regime in deriving the spin-spin interaction Hamiltonian, so ground state cooling is not necessary [61].
We can also obtain effective spin-spin interactions in equation 2 for the case when η m Ω I δ m does not hold, but in this scenario, the phonons in each mode m only return to their initial state at particular times 1/δ m . A favorable situation occurs when η m Ω I δ m for m = n, that is for all phonon modes except one. In this case, the spinspin interaction takes the form of equation 2 at times t gate = k/δ n with k an integer. This situation is usually employed when considering quantum gate operations [61]. Note that any remaining spin-motion entanglement causes decoherence in the spin state and appears as an error in the quantum simulator.
Using additional sets (two) of Raman fields we can also generate spin-spin couplings J ij ⊥ of the formĤ XX +Ĥ Y Y ∼ S x ⊗ S x + S y ⊗ S y [16,62]. For that we assume ω I is tuned close to ω ↑↓ ± ω m , a situation in which we have spin flip as well as the creation of phonons. The expression for J ij ⊥ is still given by equation 2, but with a re-defined detuning δ m = ω I − (ω ↑↓ ± ω m ). Note that the different spin-spin couplings, J ij z and J ij ⊥ , can be independently controlled by the laser parameters of the corresponding Raman fields. The corrections to the total Hamiltonian, arising from the non-commutativity of the different coupling Hamiltonians are oscillatory and can be neglected [25,63,64]. As experimental complexity increases when considering more laser beams, in Appendix D we also discuss a reduced Hamiltonian of the form H XZ =Ĥ XX +Ĥ ZZ , in which the couplingsĤ Y Y are not included. Notice that, within the perturbation theory on the gauge-invariant manifold, the strong coupling Hamiltonians of the two different cases are very similar at lowest orders. We found that, in theĤ XZ case, the underlying physics is almost unaffected when compared to the full setup. Alternatively, the simulation may be implemented by a stroboscopic sequence of "Trotter steps" in which each term in the Hamiltonian is switched on for a short time ∆t consequetively [18]. In this concatenated form, errors due to noncommutativity between the terms in the Hamiltonian scale with the discretized stepsize squared: O(∆t 2 ). This approach also has the advantage that less independent laser beams are necessary.

Single hexagonal plaquette in an N = 7 ion crystal
From equation (2) it is clear that phonon-mode spectrum is of central importance to the spatial form of the spin-spin interactions, as it determine the Lamb-Dicke parameters η j m as well as the mode detuning δ m . In single species ion crystals, in which all ions experience the same trapping frequency, the phonon modes are highly collective, which results in long range spin-spin interactions. Here, we discuss how to generate antiferromagnetic spin-spin interactions displaying a hexagonal-plaquette pattern in a 2D ion crystal. By plaquette pattern we imply that each of the six spins occupying the corners of a hexagon interact with every other spin in the same hexagon with the same strength irrespective of their inter-spin separation. First, we discuss how to generate it in a minimal setup -a crystal consisting of seven ions -and then extend to a larger crystal made of 19 ions.
The equilibrium configuration of an ion-crystal with 7 ions in a radially symmetric trap is shown in the inset of figure 2a. Although they form a triangular lattice, we picture it as a hexagonal structure with an ion in its centre. The central ion is prepared in the atomic state |1 , (see figure 1), which is immune to the Raman fields used for spin operations, but experiences an additional one-dimensional pinning optical lattice applied normal to the plane of the crystal, i.e. along the z-axis. The rest of the ions occupying the spin states form a hexagon shape. The pinning laser modifies the transversal trapping frequency of the central ion toω z = ω 2 z − ω 2 OL compared to the other ions, when the maximum of the lattice potential is focused at its equilibrium position, where ω OL is the harmonic-frequency of the local potential at the ion position due to the optical lattice. For a red detuned lattice laser and the ion trapped in an antinode, the pinning lattice relaxes the confinement of the central ion along the transversal direction and we see below how it affects the phonon spectrum. In the following, we restrict ourselves to pinning-lattice parameters such that the equilibrium positions of the ions are not altered once it is switched on.
The normal-mode spectrum consists of 14 xy (in-plane) and 7 z-(transversal) phonon modes, see figures 3a and 3b. Due to the 2D character of the crystal, the in-plane (PM) and the transverse modes (TM) are completely de-coupled. The TMs account for the oscillation of ions along the tightly confined direction, and lie on top of the energy spectrum, see figure 3a. The pinning affects mostly the TMs in two ways: , with δ 1 = 2π × 10 kHz, the detuning from the lowest TM. The slight imperfections from the BFG plaquette interactions arise from the off-resonant coupling to higher TMs. For sufficiently small value ofω z , the lowest TM mode of the 2D crystal become energetically unstable, i.e. ν 1 = 0 and leads to a structural phase transition in which the 2D character of the whole crystal has been lost [31]. For the set of parameters in figure 2 it happens wheñ ω z < ω z /2 ∼ 2π × 1.5 MHz, which restricts the frequency due to the optical lattice potential to values ω OL /ω z < 3/4. figure 2a. Note that it exhibits a hexagonal character in which all the ions in the hexagonal ring oscillate with the same amplitude and in phase. When the Raman beat frequency is tuned close to this mode, the resulting spin-spin couplings J ij z exhibit a hexagonal-plaquette pattern as shown in figure 2b. There are slight imperfections from the ideal hexagonal pattern discussed in the BFG model (equation 1) due to the offresonant couplings between the Raman fields and other TMs. These imperfections can be reduced by providing a sufficient gap between the first TM and the rest. As shown in figure 3c, this can be done by tuning the frequency of the local harmonic potential due to the pinning optical lattice at the central ion position.
If one considers a Paul trap, the anharmonic coupling between the in-plane micromotion and the transverse modes has to be taken into account, which, as recently shown [65], only introduces an overall renormalization in the ion positions without effecting the normal mode structure. While an overall shift in the transversal mode frequencies should also to be taken into account, but it has no effect on our underlying scheme MHz. The red-filled bubble indicates the lowest TM and the frequency of it is shifted to a lower value due to the presence of pinning field (shown by a black arrow). (c) The lowest transversal mode ν 1 (filled squares) and the energy difference between the lowest two transversal modes δν = ν 2 − ν 1 (filled circles) as a function ofω z are shown. Forω z < ω z /2 the mode ν 1 gets energetically unstable and the 2D nature of the crystal is lost.
for generating spin-spin interactions. Hence, we can safely neglect those effects in our calculations. In the spectrum, there exists a zero-energy in-plane mode due to the rotational symmetry of the ion crystal, corresponding to the free rotation of ion crystal in the xy plane. In real experiments, any weak stray fields will break this symmtery and the crystal will be trapped in a stable configuration.

Double hexagonal plaquette in an N = 19 ion crystal
Having established the basic idea of hexagonal plaquette interactions via an ion crystal of seven ions, we expand our discussion to two plaquettes in a larger crystal. Here, we consider a crystal of 19 ions in a radially symmetric trap. In order to create two plaquettes, we require two pinning standing wave lasers as shown in figure 1a and focus them such that the maxima of lattice potentials lie at the equilibrium positions of the two ions populated in state |1 . These two ions, as in the previous case, experiences a shallower trap along the z axis with respect to the other ions. The parameters of the pinning lattices are chosen such that the equilibrium positions are not altered. The trapping parameters are taken to be ω x = ω y = 2π × 1 MHz and ω z = 2π × 3.5 MHz and the effective harmonic frequency at the pinned ions along the z axis areω 45 MHz respectively for ion 1 and ion 2. The reason for the asymmetry in the pinning lattices parameters is discussed below. With these parameters, the two lowest TMs have the plaquette character and they comprise mostly of the oscillations of the two pinned ions. If we switch off either one of the pinning lattices, we retrieve the same physics discussed in the case for one plaquette. As an example, we switch off the pinning lattice on ion 2 and the ion 1 experiences an effective harmonic trapping frequency ofω 1 z = 2π × 2.1 MHz. The resulting lowest TM is shown in figure 4a. It accounts mostly the oscillation of ion 1 (not shown in the figure), and is well localized among the six ions surrounding it, forming a hexagonal oscillation pattern.
The resulting spin-spin interactions, when the Raman fields are detuned close to this mode, exhibit the hexagonal-plaquette pattern (see figure 4b) similar to what we have discussed in section 2.2 for the case of 7-ion crystal. This point constitutes one of the main results of our paper: even in a larger crystal, by engineering phonon modes we can create localized spin-spin interactions, which in our case exhibit a hexagonal pattern. Compared to the smaller crystal, here the imperfections are slightly higher. This can be understood from the structure of the plaquette eigenmode itself: the amplitude of oscillation is not the same for all the six ions in the hexagon (see figure 4a), which is contrary to the N = 7 case. This arises due to the asymmetry in the number of nearestneighbor ions for each ion in the hexagonal ring. A second reason is that the number of modes increases linearly with the system size, but in our case the coupling to those modes are well suppressed by the large detunings δ m from the corresponding m th mode. Note that if we pin ion 2 instead of ion 1, we create a hexagonal plaquette around the second ion.
Two plaquettes case: To create two plaquettes, we switch on both pinning lattices. In this case, the two lowest phonon modes are used to engineer the plaquette interactions. If we have identical pinning lattices, the two modes become degenerate. Then, dressing them with a single pair of Raman fields may generate an arbitrary superposition of phonon states: c 1 b 1 ± c 2 b 2 , where b 1,2 are the eigen vectors of the two lowest TMs and |c 1,2 | 2 provide us the population of the respective phonon modes. The resulting spin-spin interactions do not possess the plaquette character. Hence, we introduce an asymmetry between the parameters of the pinning lattices such that the degeneracy is lifted as well as the modes are well separated in the spectrum. Then, we use two different pairs of Raman fields to induce spin-spin couplings, in which one is detuned near the mode b 1 and the second one is near the mode b 2 with the same detuning δ. Locally, we can control the interaction strengths in each plaquette by tuning the effective coupling strengths due to two Raman fields, i.e., Ω 1 I and Ω 2 I ; in the particular example shown in figure 5b we choose Ω 2 I /Ω 1 I = 1.1 to obtain the same interaction strengths in both plaquettes. The equilibrium configuration of the 2D ion crystal and the eigenvectors of The hopping parameter J ⊥ is generated by additional Raman fields, again mediated through transversal modes. The nature of J ⊥ is also vital in determining the ground state properties as well as the magnetization dynamics of the emulating spin Hamiltonian. As short-range hopping is preferable, e.g. of nearest-neighbor type [8], we identify three lowest TMs in the spectrum, labelled as m =1, 2 and 3, see figure 6a. They lie just above the plaquette modes and are also well separated from them, guaranteeing that the Ising part is not influenced when addressing with the Raman fields. Note that these modes are shared among the two corner-sharing triangles at the centre of the crystal, and hence the resulting hopping dynamics may introduce very interesting interplaquette spin dynamics, for example ring exchange among the four ions for Ω 2 I /Ω 1 I = 1.1 and δ = 2π × 20kHz. Note that ion 6 is shared by both the plaquettes.

Many-body Physics and Quantum Magnetism in a double-plaquette quantum simulator
In this section we discuss the many-body physics associated with the above mentioned spin systems with hexagonal-plaquette interactions. In particular, we focus on the ground state properties and magnetization dynamics. First, we consider the plaquette interactions alone, and address the role of the imperfections on the interaction patterns at the classical level. Then, we show how gauge-invariance, embodied by the magnetization conservation on each plaquette, is affected by the introduction of quantum fluctuations, and compare the realistic scenario with the ideal BFG model. Finally, we show how interesting many-body dynamics can be observed in a minimal system involving two plaquettes, within reach of state-of-the-art experiments with trapped ion crystals, by exploring the dynamics of a single charge on top of the gauge-invariant background.

Ground state properties
Single plaquette case: In the classical limit (J ⊥ = 0) of the HamiltonianĤ (equation 1) in a single plaquette, the ground state manifold is featured by 20-fold degenerate states. Each of them obeys the constraint i∈ S i z = 0, i.e., the total magnetization around the hexagon is zero or, equivalently, three spins are "up" and the other three spins are "down". All such configurations are shown in figure 7. The constraint is akin to the wellknown ice rule "2-in 2-out" for the orientation of four dipoles arranged at the corners of a tetrahedron in a crystal ice. For the spin-version of the ice setup, so called the "spinice", the ice rule reads as "2-up 2-down" [3,43] and can be mapped into a square lattice in which around each vertex four spins are arranged according to ice-rule. There are six such possible spin-configurations. This six-vertex or ice model required that all the four spins around the vertex interact with the same strength, similar to the BFG model in which all the spins around the hexagon have the same pairwise interaction strength.
In our ion-setup, there are imperfections from a perfect plaquette pattern as we have discussed above, which lift the 20-fold degeneracy partially. However, as imperfections are mostly involving spins within or on adjacent plaquettes, the investigation of a two plaquette system is required to access their effects.  Double plaquette case: The number of degenerate ground states increases with the system size, and the degree of degeneracy is 200 for two plaquettes with the Ising part of the BFG model. In the ion setup, as in the previous case, the imperfections lift this degeneracy partially, but the energey levels are closely spaced in the energy spectrum as is evident from figure 8a. In our particular example for the ion setup, the classical ground states are featured by four degenerate sates and even a very low temperature can lead to a classical order-by-disorder phenomenon in which the system collapse into one of the ground state configuration. To show the contrast, we compare the above results with a nearest-neighbour Ising model in two-plaquettes, with no inter-plaquette spinspin interactions, the ground states are doubly degenerate (antiferromagnetic states) and any excited state, resulting from a single spin flip, leads to an energy cost of J z , see figure  8a. This implies that our ion setup for the implementation of plaquette interactions is indeed a promising candidate for emulating frustrated quantum magnetism and en route to quantum spin liquid once scaled up to larger systems.
In the presence of quantum dynamics the global magnetization, M = iŜ i z is preserved, but the magnetization per palquette is not. We define the following operator to quantify the magnetization per plaquette, where N is the number of hexagonal plaquettes (= 2 in our case). The magnitude of G , the expectation value taken in the ground state of the spin-Hamiltonian, gives us the measure of the admixture from states outside the classical ground state manifold obeying the plaquette constraint. In practice, this determines the quality of gauge-invariance as a function of quantum fluctuations: indeed,Ĝ can be viewed as an effective generator of a local (gauge) symmetry, which commutes with the effective Hamiltonian in the strong coupling J z J ⊥ limit [10]. With no hopping (J ⊥ = 0),Ĝ = 0, and it increases as J ⊥ increases as shown in figure  8b. The small values of Ĝ even for sufficiently large values of J ⊥ shows the stability of the low energy manifold with a plaquette constraint against quantum fluctuations. In both cases, as expected from perturbation theory arguments, Ĝ ∼ (J ⊥ /J z ) 2 close to the classical limit. We note that the only effect of the imperfections is quantitative: some gauge variant states not satisfying the plaquette constraint pay a smaller energy penalty compare to the ideal BFG case. However, their effect is relatively small even for intermediate coupling strengths. The small values of Ĝ even for a sufficiently large values of J ⊥ guarantees us the the stability of gauge-invariant low energy manifold. The upper axis for the Rabi frequency (corresponding to the J ⊥ in the lower axis) is estimated for the particular case when Ω I = 2π × 500 kHz, δ m = 2π × 20 kHz and δ ⊥ = −2π × 60 kHz as discussed in section 4. In calculating Ĝ for the BFG model, we restrict the hopping in the five spins located at the intersection of the two plaquettes, which makes sense when considering the system in the thermodynamic limit and also for the hopping pattern shown in figure 6d for the real situation.

Plaquette-magnetization dynamics: deconfined charge excitation
In this section we examine and compare the dynamics of magnetization in each plaquette for an ion setup with plaquette interactions with the BFG model, for different initially prepared spin configurations. The local plaquette magnetizations are defined as M 1 = i∈ 1 S i z and M 2 = i∈ 2 S i z respectively for plaquettes 1 and 2. To address the stability of gauge constraint against spin hopping, we study the time propagation starting from one of the classical degenerate ground state, which obeys 3-in 3-out rule. As expected, the hopping dynamics hardly affected the local plaquette magnetizations even for large hopping parameter J ⊥ . A particular example with J ⊥ /J z = 0.2 is shown in figures 9a and 9b for the ion setup and the BFG model, respectively, and they are found to exhibit the same behaviour. The dynamics become more interesting when we introduce a charge in one of the plaquttes by flipping a single spin in the respective plaquette. As an example we consider an initial state in which M 1 = 0 and M 2 = −1 and the inter-plaquette spin hopping couples this state to M 1 = −1 and M 2 = 0 and hence the state oscillates between them. This has a clear interpretation in gauge theoretical language: the "doped" plaquette represents a charge moving on top of a gauge-invariant background. This scenario can be compared to a case where a pair of quasi-particle and quasi-hole (fractionalized charges or excitations) is created by a single spin flip in a square lattice, which supports an exact incompressible quantum liquid. The particlehole pair doesn't form a bound pair but is in a deconfined phase [66]. The same is also predicted for a classical Ising model in pyrochlore lattice in which the macroscopic ground state degeneracy leads to the deconfinement of monopole excitations [67]. This is one of the intriguing point of this paper from a frustrated magnetism picture, arising with a question on the stability of this deconfined phase, interestingly, which could be qualitatively addressed with a small size ion-crystal quantum simulator as we show here.
The many-body physics associated with the reduced Hamiltonian H XZ mentioned in section 2.1 is discussed in Appendix D and found to possesses the same properties as that of the full Hamiltonian.
4. Implementation with a 2D crystal of 40 Ca + As an example, we consider a crystal of 40 Ca + ions. Since this ion does not posses a nuclear spin, we employ the Zeeman and metastable electronic states, rather than hyperfine ground states. We encode the spins in electronic states of the ion, | ↓ = |S 1/2 , m j = 1/2 and | ↑ = |D 5/2 , m j = 3/2 [68]. Inititialization can be performed by optical pumping followed by addressed coherent laser pulses that prepare each ion in a designated spin state. In this way the ions that will be pinned can also be prepared in the state |1 = |S 1/2 , m j = −1/2 . The pinning can be done using e.g. a λ pin = 532 nm retro-reflected laser with a power of ∼ 50 mW, focussed to a waist of 2 µm. For an estimate, in which we are neglecting the polarizabilities of all transitions except the dominant S → P transition around 397 nm, this results in an anti-trapping potential of trap-frequencies (0.08,0.08,1.3) 2π MHz in a node of the standing wave. With a Paul trap of transverse frequency ω ⊥ = 2π 3 MHz, this results in a local reduction of the trapfrequency of ∼ 10 % to 2π 2.7 MHz as in the example in section 2.2. The trap-frequency in the other two directions is hardly affected (it is reduced by ∼ 2π 3 kHz). The lines with bubbles are the case for which the initial state is such that one of the plaquette carries a charge (in this particular example a negative charge in the second plaquette, i.e., M 2 = −1) and the other plaquette has zero charge (or zero magnetization, M 1 = 0). In time the charge oscillates between the two plaquettes indicating a de-confined phase.
Errors in the realization of the BFG model arise from spontaneous scattering from the pinning laser, wich is reduced because of the large detuning and because the ion is placed in a node of a wave. We expect a scattering rate of ≤ 1 s −1 as an upper bound. The ions are separated by about 5 µm, which is much larger than the waist of the pinning laser, such that the ions in the plaquette are not affected by it.
From an experimental point of view, it is convenient to perform a change of basis to exchange σ z and σ x , as spin-spin interactions in the latter require less laser power [69], whereas we are interested in the regime J z > J ⊥ . Spin-spin interaction terms involving σ y can be ommitted in a first implementation and this results in non-resonant S + ⊗ S + and S − ⊗ S − terms that are suppressed by the Gauss law, shown in Appendix D. The J z plaquette spin-spin interaction can now be induced by a Mølmer-Sørensen scheme [61], operating on the optical qubit at 729 nm [68,70]. In particular, a bichromatic light field with components around ω ↑↓ ± ω m ± δ m induces S x ⊗ S x interactions [69]. Rabi frequencies of Ω I = 2π 500 kHz and a detuning δ m = 2π 20 kHz result in (plaquette) spin-spin interactions of J z 2π 18.5 kHz. The pinned ion would have a transient maximal amplitude of oscillation of around 20 nm λ pin , such that the harmonic approximation is justified for the pinning laser. The fidelity of the gate for σ z is calculated semi-classically and its lower bound is about 99.7%, limited by remaining spin-motion entanglement in far detuned phonon modes. For details regarding the fidelity calculation see Appendix C. To realise J ⊥ /J z 0.2, we require a Rabi frequency for the hopping spin-spin interactions of about Ω I 2π 20 kHz with δ ⊥ = −2π 60 kHz as discussed above.
For ions with a nuclear spin, such as 9 Be + or 171 Yb + , the three required long-lived states can be encoded in the hyperfine ground states and the spin-spin interactions can be induced by Raman lasers, which can reach comparable coupling strengths as for the electronic state encoding descibed above [71]. For Yb + strong optical pinning has been shown recently [55].
The readout of the quantum simulation run is performed by detecting statedependent laser induced resonance fluorescence, where the state | ↓ will appear bright while sites with ions in | ↑ remain dark as the state does not scatter photons. Detection fidelities exceed 99% in ion trap experiments routinely. Note that the detection can easily resolve individual ions, such that the outcome can be immaged on a ccd chip in one picture, where populations, and correlations may be readily revealed. This allows not only for obtaining the global magnetization but also to determine correlations of spins such that the gauge-invariant dynamics can be observed in the plaquettes.

Conclusion and Outlook
In conclusion, we have shown that using phonon mode shaping it is possible to simulate exotic spin-spin interactions in an ion crystal which have far reaching consequences in the context of strongly correlated systems, in particular frustrated magnetic systems in which low energy manifold is described by an emergent dynamical gauge field. The mode-shaping is accomplished by (anti-) pinning standing wave light fields, which modify locally the ion trap and affect its motional dynamics. Alternatively, we can selectively excite ions to a Rydberg state [57,58] in which those ions experience a different trapping potential compared to the ground state ions due to state-dependent atomic polarizabilities, see Appendix E. Additional microwave fields could be employed for fine-tuning [59,60] the trapping field of the ions occupying the Rydberg state. This leads to a design of the tranversal mode structure in the planar crystal, similar to that with optical pinning forces.
From a quantum magnetism point of view, we have discussed the ground state and magnetization dynamics of the corresponding spin Hamiltonian of the proposed trapped ion implementation and the results have been compared with that of the idealized BFG model, within a small setup of one and two plaquettes. The results show the excellent agreement between the physics in an ion based quantum simulator and that of the BFG model. The deconfined dynamics of a single charge excitation in a two plaquette system, arising from the ring exchange, is akin to the deconfined monopole excitations in a 3D pyrochlore spin-ice, arising from the large macroscopic ground state degeneracy. In general, the current studies would open up a completely new aspect of studying spin-Hamiltonians using ultra cold ion crystals, in particular once scaled up, will be the prime quantum simulator for emulating topological quantum spin liquids or resonatingvalence-bond states. quantum gates based on an effective laser field being interacting with many ions. The electric field from the laser beams gives rise to a Stark shift for each internal states, and a corresponding electric-dipole force on each ion. These state dependent forces are typically used in trapped-ion quantum computation [61,74]. The time evolution of the system is then obtained using the Magnus expansion [75] for the time evolution operator where H(t) is the time-dependent Hamiltonian. In our case, with sufficiently large detuning from the motional sideband [16] and assuming the same phase for all ions, the long-term time evolution is dominated by a term linear in time, t and is of the form This realizes an Ising spin-Hamiltonian of the formĤ z = i,j J ij z S i z ⊗ S j z with spin-spin couplings J ij z between i th and j th ions. Note that the spin-1/2 operators (S α ) defined in equation 1 are related to the Pauli spin-1/2 matrices (σ α ) through S α = σ α /2.

Appendix B. Ion crystal: Equilibrium positions and vibrational spectrum
In a 3D setup, the potential experienced by the ions is The first term is the harmonic confinement due to external fields and the second term accounts for the Coulomb interaction between the positively charged ions. We choose the length scale as l 3 x = (Z 2 e 2 )/(4π 0 M ω 2 x ) and in the dimensionless form the potential reads as where λ z = ω z /ω x , λ y = ω y /ω x . In the following we absorb l −1 x in the position coordinates such that x i , y i and z i are dimensionless. The equilibrium positions are calculated by solving the equations for each ions. This leads to solving the coupled algebraic equations of the form where α = {x, y, z} and R 0 (i, j) = (x i −x j ) 2 +(y i −y j ) 2 +(z i −z j ) 2 . Once the equilibrium positions are obtained, we can calculate the eigen vectors b m and eigen values ν m of the phonon modes by the exact diagonalization of the Hessian matrix constructed out of the following second order derivatives: with α, β ∈ {x, y, z} and α = β. after we left out the zero point oscillation energy. This relation would provide us a rough estimation for the number of phonons in the system for a particular mode m. In addition, the overlap between the initial phonon state and the maximally displaced state would give us the fidelity loss arising from this particular mode. We require this loss would be well below 1%. Taking the Gaussian ground state solution of the harmonic oscillator, we get the overlap as a function of the displacement z 0 as O m = exp(−z 2 0 /8q 2 m0 ). (C.5) The fidelity loss L m = 1 − O m due to a given mode m, as a function of the amplitude of oscillation is shown in Fig. C1, and one can see that the larger the displacements, the lesser the fidelity is. In order to keep the fidelity loss with in 1%, the displacement z/q m0 required to be ≤0. 30. The net fidelity loss is then obtained by the summing over that of all modes. Appendix E. Rydberg excitations in an ion crystal: atomic polarizibility and modified trap frequency Instead of using the pinning lattices on ions in state |1 (see main text), we can as well excite or weakly dress those ions to a Rydberg state. Due to the state dependent atomic polarizabilities, the Rydberg excited ions experience different trapping frequencies compared to those occupying the low-lying atomic states. This modifies the phonon spectrum of the ion crystal, identical to the situation in which the pinning lattices are present. Recently, it has been shown [59,60] that ions excited to the nP -Rydberg state experience an additional radial potential, V a (r) ≈ −e 2 α 2 P nP r 2 , where the polarisability P nP ≈ −0.25 × n 7 atomic units, with n, the principal quantum number of the Rydberg state. Effectively, an ion in this Rydberg state experiences a tighter confinement compared to an ion occupying a low-lying state, but with additional microwave fields, by coupling to a nearby Rydberg state, say n S, we can freely tune the trapping frequency of the Rydberg excited ion. It can be increased, decreased or even made equal to that of ions occupying the qubit states | ↑ and | ↓ . The dressed states |± = N ± (C ± |nP + |n S ) polarizabilities read as P ± = N 2 ± C 2 ± P np + P n S (E.1) with P n S > 0, the parameter C ± depends on the microwave parameters such as detunings and Rabi frequencies and N ± are the normalization constants. More details on the microwave-Rydberg approach can be found in [60]. P ± are easily controlled via microwave field paremeters, and hence the local trapping of Rydberg ions. The only limit in this approach is set by the life time of Rydberg states for gate operations and it can be augmented by going either to a higher Rydberg level or by weakly dressing to the Rydberg state.