Quantum Magnetism of Spin-Ladder Compounds with Trapped-Ion Crystals

The quest for experimental platforms that allow for the exploration, and even control, of the interplay of low dimensionality and frustration is a fundamental challenge in several fields of quantum many-body physics, such as quantum magnetism. Here, we propose the use of cold crystals of trapped ions to study a variety of frustrated quantum spin ladders. By optimizing the trap geometry, we show how to tailor the low dimensionality of the models by changing the number of legs of the ladders. Combined with a method for selectively hiding of ions provided by laser addressing, it becomes possible to synthesize stripes of both triangular and Kagome lattices. Besides, the degree of frustration of the phonon-mediated spin interactions can be controlled by shaping the trap frequencies. We support our theoretical considerations by initial experiments with planar ion crystals, where a high and tunable anisotropy of the radial trap frequencies is demonstrated. We take into account an extensive list of possible error sources under typical experimental conditions, and describe explicit regimes that guarantee the validity of our scheme.

The transition from single-particle to many-body quantum systems yields one of the most interesting concepts of physics: emergence. As emphasized by P. W. Anderson [1], the laws that describe the collective behavior of interacting many-body systems can be fundamentally different from those that govern each of the individual particles. This effect leads to some of the most exotic phenomena of condensed-matter physics in the last decades [2]. Unfortunately, the high complexity of many-body systems turns our endeavor to understand such emergent phenomena into a fundamental challenge for both experimental and theoretical physics. Since exact analytical solutions seldom exist, even for oversimplified models, one is urged to develop efficient numerical methods. An alternative to this approach are the so-called quantum simulations [3], which make use of a particular quantum system that can be experimentally controlled to a high extent, in order to unveil the properties of a complicated interacting many-body model.
There are two different strategies for the quantum simulation of many-body systems [4], the so-called analog and digital quantum simulators (QSs). On the one hand, an ideal analog QS would be a dedicated experimental device where one can prepare the quantum state, engineer a Hamiltonian of interest, and measure its properties. On the other hand, a digital QS aims at reproducing the dynamics of any given Hamiltonian by concatenating a set of available quantum gates. Independently of their particular experimental implementation, these types of QSs would be capable of exploring models from very different areas of physics, ranging from condensedmatter to high-energy physics. Let us remark that this enterprise benefits directly from the development of architectures for quantum-information processing [5]. In this work, we focus on one of these architectures: laser-cooled atomic ions confined in radio-frequency traps [6]. Here, some remarkable quantum simulations have already been accomplished in experiments, either in the digital [7] or analog [8][9][10][11][12] versions. In particular, we shall concentrate on analog QSs, where a arXiv:1205.0341v2 [quant-ph] 26 Sep 2012 number of theoretical proposals already exist [13]. Some of these proposals target the quantum simulation of relativistic effects [14], quantum spin models [15][16][17][18], many-body boson systems [19], spin-boson models [20], or theories of quantum transport and friction [21].
A direction of research that is being actively explored is the QS of magnetism [15], which yields a playground for a variety of cooperative quantum many-body effects. Even though the phenomenon of magnetic ordering was already known to the ancient greeks, the understanding of its microscopic origin could only be achieved after the development of quantum mechanics. A particular topic of recurring interest is the fate of the magnetically-ordered phases in the presence of quantum fluctuations [22], which get more pronounced as the dimensionality of the system is reduced. Such quantum fluctuations may destroy the long-range order, favoring paramagnetic phases and triggering the so-called quantum phase transitions [23]. Alternatively, they may stabilize more interesting phases, such as the long-sought two-dimensional quantum spin liquids, which have connections with the theory of high-temperature superconductivity [24]. The advent of a QS for quantum magnetism would allow for an unprecedented experimental realization of these phases, addressing puzzles that remain unsolved due to their great complexity.
A representative model that exemplifies the usefulness of analog QS is the one-dimensional quantum Ising model [25], which describes a chain of interacting spins subjected to a transverse magnetic field. This model has traditionally been considered as a cornerstone in the theory of quantum phase transitions [23]. However, the vast majority of low-dimensional materials realize instances of the so-called Heisenberg model (see e.g. [26]). In fact, the demanding requirements to explore the Ising magnetism in real materials (e.g. precise one-dimensionality, strong Ising anisotropy, and weak exchange couplings matching the available magnetic fields) have postponed its observation until the recent experiments with CoNb 2 O 6 [27][28][29][30]. We stress that prior to this experiment, trapped-ion QSs had already targeted the onset of Ising criticality [8,10]. More recently, Ising interactions in a large collection of ions have been observed in Penning traps [12], which in combination with the recent results for neutral atoms [31], yield the unique possibility of tailoring the microscopic properties of the quantum Ising magnet (e.g. couplings, dimensionality, geometry).
In this manuscript, we investigate the capabilities of trapped ions as QSs of frustrated quantum Ising models (FQIMs). Here, frustration arises from the impossibility of minimizing simultaneously a set of competing commuting interactions (i.e. classical frustration) [32], whereas the quantum fluctuations are introduced by a transverse magnetic field. We remark that, with the exceptions of three-dimensional spin ice [33] and some disordered spin-glass compounds [34], most of the frustrated materials correspond once more to Heisenberg magnets [35]. This would turn the extensive research on lowdimensional frustrated Ising models [36,37] into a pure theoretical enterprise. Fortunately, the seminal experiment [9] has proved the contrary by demonstrating that the physics of small frustrated networks [17] can be explored in trapped-ion labo-ratories. In a recent work [18], we have studied a different approach to control the degree of frustration in ion crystals, which works independently of the number of ions and is thus amenable of being scaled to large systems. In this manuscript, we elaborate on that proposal to develop a versatile quantum simulator for a variety of spin ladder geometries.
A quantum spin ladder is an array of coupled quantum spin chains. Each of these spin chains is usually known as a leg of the ladder, whereas the couplings between them define the ladder rungs. Let us note that the study of antiferromagnetic Heisenberg ladders has a long tradition inspired by the experiments with insulating cuprate compounds [38], which become high-temperature superconductors after doping [39]. These strongly-correlated systems lead to fascinating and thoroughly-studied phenomena [40]. In contrast, the quantum Ising ladders have remained largely unexplored, probably due to the absence of materials that realize them. Our work discusses a proposal to fill in such gap, which is based on techniques of current trapped-ion experiments.
This article is organized as follows. In Sec. II, we describe the properties of the ladder compounds formed by a collection of ions confined in a radio-frequency trap. We discuss the experimental conditions leading to a particular vibrational band structure, which shall be exploited to explore the physics of magnetic frustration. In section III, we support the above conclusions by the numerical simulation of a particular case: the trapped-ion zigzag ladder. A discussion of the experimental feasibility of our scheme with state-of-the-art setups is also presented in Sec. IV. Here, we also present initial experiments regarding the appropriate trap design for the QS, and discuss possible sources of error. In Sec. V, we list the many-body phenomena that can be addressed with the proposed QS, ranging from ordered phases of quantum dimer models, to disordered quantum spin liquids. In Sec. VI, we address in detail the phase diagram of the dipolar J 1 -J 2 quantum Ising model. Finally, we present our conclusions in Sec. VII. In the Appendixes A, B and C, we discuss technical details regarding both the spin-dependent dipole forces and spontaneous decay, the micromotion, and the thermal secular motion. We describe conditions under which they do not affect the proposed QS.

A. Geometry of the trapped-ion ladders
An ensemble of N atomic ions of mass m, and charge e, can be trapped in a microscopic region of space by means of radio-frequency fields [6]. This system is described by where {ω α } α=x,y,z represent the effective trapping frequencies (see also Sec. IV and Appendix B), and {r i , p i } N i=1 are the positions and momenta of the ions. Note that we use gaussian units andh = 1 throughout this text.
With the advent of laser cooling, it has become possible to reduce the temperature of the ions to such an extent that they crystallize. This gives rise to one of the forms of condensed matter with the lowest attained density, ranging from clusters composed by two ions [41], to crystals of several thousands [42]. Interestingly, these crystals can be assembled sequentially by increasing the number of trapped ions one by one, so that the aforementioned transition to the many-body regime [1] can be studied in the laboratory. Besides, the geometry of the crystals can be controlled experimentally by tuning the anisotropy of the trapping frequencies [43], which opens a vast amount of possibilities for the QS of magnetism.
The ion equilibrium positions {r 0 i } N i=1 are determined by the balance of the trapping forces and the Coulomb repulsion. Formally, they are given by ∇ ∇ ∇ r i V = 0, where V contains the trapping and Coulomb potentials. By introducing the unit of length l z = (e 2 /mω 2 z ) 1/3 , the dimensionless equilibrium positionsr 0 i = r 0 i /l z follow from the solution of the system where we have introduced the anisotropy parameters κ α = (ω z /ω α ) 2 . As shown in the experiments [44], by increasing the parameters κ x = κ y towards a maximum value κ max = 1, the geometry of the ion crystal undergoes a series of phase transitions starting from a linear chain, via a two-dimensional zigzag ladder, to a three-dimensional helix. We remark that these structural phase transitions have been the subject of considerable interest on their own [45,46]. We have recently investigated the possibility of setting a large anisotropy between the radial trap frequencies 1 ≥ κ x κ y , such that the ion crystals get pinned to the xz-plane [18]. This requires a modification of the symmetric electrode configuration of the usual linear traps (see the details in Sec. IV). Besides, the parameter κ x can be increased by a bias voltage, yielding a variety of geometries. For instance, above a certain value κ c,2 , the linear chain transforms onto a 2-leg ladder (see Figs. 1(a)-(b)). By solving numerically the system of equations (2), we observe that the 2-leg zigzag ladder first evolves into a 3-leg ladder as κ x is increased, then yields a 4-leg ladder [see Figs. 1(c) and 1(f)], and so on. Hence, there should be a sequence of structural transitions at the critical values κ c,2 < κ c,3 < κ c,4 < · · · < κ c,n l , where n l determines the number of legs of the trapped-ion ladders. The rungs of these ladders correspond to the diagonal links, yielding a collection of bond-sharing triangular plaquettes. In Sec. IV, we describe a method to build ladders from corner-sharing triangles (see Figs. 1(d)-(e)), which widens the applicability of our QS.
Eventually, when κ x → 1, the crystal must correspond to a two-dimensional bond-sharing triangular lattice, or cornersharing Kagome lattice, with ellipsoidal boundaries. Accordingly, not only can we control the number of legs of the trapped-ion ladder, but also explore the crossover from quasione-dimensional physics to the two-dimensional realm. With respect to the FQIMs, this two-dimensional limit is very appealing due to its connection to quantum spin liquids [24], and quantum dimer models [37]. However, as happens for Heisen-  Figure 1. Trapped-ion ladders: Self-assembled geometries for N trapped ions as a function of the anisotropic trapping frequencies ω y ω x , ω z : (a) Linear chain for N = 50 and κ x = 4 · 10 −4 , (b) Twoleg zigzag ladder for N = 51 and κ x = 1 · 10 −2 . (c) Three-leg ladder of bond-sharing triangles for N = 52 and κ x = 4 · 10 −2 , (d) Two-leg ladder of corner-sharing triangles for N = 52. From this geometry, it is possible to obtain the two following ladders by selectively hiding some of the ions (grey crosses) following the methods outlined in Sec. IV. (e) Three-leg ladder corresponding to a Kagome stripe of corner-sharing triangles and hexagons for N = 52, (f) Four-leg ladder of corner-sharing triangles for N = 53 and κ x = 8 · 10 −2 .
berg magnets [38], the ladder compounds may already contain fascinating phenomenology. Let us remark that the precise knowledge of the crystal plane is a fundamental advantage for the quantum simulation of FQIM [18]. Prior to the discussion of the effective FQIM, we will describe the collective vibra-tional modes in these ladder compounds, since they shall act as mediators of the magnetic Ising-type interaction.
At this point, it is worth commenting on the possibility of engineering the lattice positions of the ion crystal by means of micro-fabricated electrode arrangements, (e.g. surface traps as proposed in [47]). If the technical problems related to the anomalous heating observed close to the electrodes are overcome, these new generation of traps could be combined with the proposed QS to study frustrated spin models in arbitrary lattices without the complication of micromotion. However, to keep within reach of the current technology, we focus on the more standard Paul traps where experiments on quantum magnetism have already been performed [8][9][10].

B. Collective vibrational modes of the ion ladder
Due to the Coulomb interaction, the vibrations of the ions around the equilibrium positions, r i = r 0 i + ∆r i , become coupled. By expanding the Hamiltonian (1) to second order in the displacements ∆r i , one obtains a model of coupled harmonic oscillators corresponding to the harmonic approximation In this expression, the coupling matrix between the harmonic oscillators can be expressed in terms of the mutual ion distance r 0 i j = r 0 i − r 0 j , and the quadrupole moment for each pair of ions Q For the particular ladder geometries in Fig. 1, we identify two types of vibrational excitations, the so-called transverse modes that correspond to the vibrations of the ions perpendicular to the ladder (i.e. y-axis), and the planar modes that account for the coupled vibrations within the ladder (i.e. x, z axes). The transverse modes are decoupled from the planar vibrations, and are described by a set of coupled oscillators . This term amounts to the Einstein model of individual lattice vibrations, where the renormalization of the frequencies is caused by the mean-field-type interaction of one ion with the rest of the ion ensemble. In our case, we must also consider the coupling between distant oscillators, whose magnitude relative to the trapping frequencies scales as (e 2 /l 3 z )/(mω 2 y ) = κ y 1. Hence, the transverse vibrations are described by a set of harmonic oscillators with weak couplings that decay with a dipolar law.
The planar vibrations are more complex since the motion along both axes, x and z, becomes coupled through the Coulomb interaction. In this case, the Hamiltonian is il is expressed in terms of the dimensionless equilibrium positions. Since we have assumed that κ x κ y , the planar vibrations are coupled more strongly than the transverse ones. In the last term of this Eq. (6), we observe how the nondiagonal terms of the quadrupole, Q xz i j = 0, are responsible for the coupling of the motion along the x and z axes of the ladder.
So far, we have reduced the original Hamiltonian (1) to a pair of quadratic boson models (5)- (6). These can be exactly solved by introducing the canonical transformations where a † n , a n stand for the N creation-annihilation operators for the quantized excitations (i.e. phonons) of the transverse modes, and b † n , b n are the the corresponding 2N operators for the planar modes. In these expressions, M ⊥ i,n (M i,n ) determines the amplitude of the transverse (planar) oscillations of the ion at site i due to the collective vibrational mode labeled by n, such that Ω ⊥ n (Ω n ) are the corresponding normal-mode frequencies. Formally [48], M ⊥ i,n (M i,n ) are given by the orthogonal matrices that diagonalize the second order Coulomb couplings in Eqs. (5)-(6), or equivalently are dimensionless couplings. These equations (8) must be solved numerically with the previous knowledge of the equilibrium positions (2), and yield the following quadratic phonon Hamiltonian n a † n a n + where the index n labels the normal modes with an increasing vibrational frequency Ω ⊥ n+1 > Ω ⊥ n , Ω n+1 > Ω n . a n N 2N Vibrational modes of the ladder compounds: (a) Schematic representation of the phonon branches, whereby the transverse vibrational frequencies Ω ⊥ n span around the trap frequency ω y , and are situated far away from the planar phonon branch Ω n , which spans around ω x , ω z . (b) Laser beam arrangement (inset) for a spindependent dipole force with an effective wavevector k L in the xyplane slightly tilted from the y-axis. The effective frequency ω L is tuned above the resonance of the transverse modes (see (a)), such that the detunings fulfill δ δ ⊥ .
We now discuss a qualitative picture for the phonon branches valid for all the different ladders. As argued above, the transverse phonon modes are weakly coupled due to the small parameter κ y 1. By inspecting Eq. (8), one realizes that this condition will lead to a transverse phonon branch with a small width around the trap frequency ω y (see Fig. 2(a)). This property is not fulfilled by the planar modes Ω n ∈ [Ω 1 , Ω 2N ], which are more strongly coupled κ x κ y , and will thus present a wider branch around ω x , ω z . However, due to the frequency anisotropy ω y ω x , ω z , the planar phonon branch will always be situated far away from the transverse-mode frequencies, even if its width is considerably larger. This property will turn out to be essential for the QS of frustrated spin models [18], whereby the role of the spin is played by two electronic levels of the ion, and the interactions are mediated by the collective vibrational excitations.
The main idea is that such a clustered phonon spectrum will allow us to use a pair of laser beams that only couple the spins to the transverse phonon modes, even when the radiation resulting from the interference also propagates along the plane of the ladder ( Fig. 2(b)). The interest of this idea is two-fold. On the one hand, the transverse phonon modes are the ideal mediators of the spin-spin interactions due to their higher insensitivity to ion-heating and their lower contributions to the thermal noise in the QS, as compared to the planar modes. On the other hand, by using a laser configuration with a component along the rungs of the ladder, we can exploit the ratio of its effective wavelength with the ion mutual distances in order to tailor the sign (ferromagnetic/antiferromagnetic) and the magnitude of the spin-spin interactions anisotropically (i.e. depending on the direction joining the pair of ions). As shown below, this opens the possibility of realizing a versatile QS of frustrated quantum spin ladders, which is amenable of being scaled to larger ion numbers and a variety of geometries.

C. State-dependent dipole forces
The use of the collective vibrational modes as a common data bus to perform two-qubit gates can also be understood in terms of phonon-mediated spin-spin interactions [49], which opens the route towards the QS of quantum magnetism [15,16]. We now introduce the key ingredient for such QSs, namely, a spin-phonon coupling that originates from a laser-induced state-dependent dipole force [50,51].
So far, our discussion applies to all ion species and radiofrequency traps. However, in order to exploit the phonons as mediators of a spin-spin interaction that is sufficiently strong, the trap frequencies for each particular ion species must be tuned such that l z = (e 2 /mω 2 z ) 1/3 ≈ 1-10 µm, and ω y ω x ≥ ω z . This typically constraints the order of magnitude of the trap frequencies to the range ω z /2π, ω x /2π ≈ 0.1-1 MHz, and ω y /2π ≈ 1-10 MHz. Moreover, since we aim at a flexible control of the anisotropy of the spin couplings, we shall focus on the ion species that allow for a two-photon lambda scheme to implement the spin-phonon coupling. Hence, our discussion will be specific to singly ionized alkaline-earth ions, either with a hyperfine structure 9 Be + , 25 Mg + , and 43 Ca + , or with a pair of Zeeman-split levels such as 40 Ca + and 24 Mg + [52].
In Fig. 3, we represent schematically the atomic energy levels of such ions, which have a single valence electron in the orbital n 2 S 1/2 that can be optically excited to n 2 P 1/2 , n 2 P 3/2 via a dipole-allowed transition. Here, we use the standard notation n 2S+1 L J , with n as the principal quantum number, and S, L, J as the spin, orbital, and total electronic angular momentum. Depending on the nuclear spin I of the particular ion, the ground-state manifold will be split into a pair of Zeeman states (I = 0), or a set of hyperfine levels (I = 0). We select two of such states |↑ i , |↓ i , which are separated by an energy gap ω 0 , to form the effective spins of our QS (see the inset of Fig. 3). For hyperfine spins, the energy splitting is ω 0 /2π ≈ 1-10 GHz, whereas for Zeeman spins its order of magnitude depends on the external magnetic field ω 0 /B 0 ≈ 2π×10 GHz/T. Accordingly, the resonance frequencies usually lie in the radio-frequency/microwave regime. However, it is not customary to use radio-frequency/microwave radiation to couple the spins to the collective vibrational modes, since its wavelength is much larger than the typical ion oscillations (but see the recent experimental progress [53]). A possible alternative is to use a pair of laser beams with optical frequencies ω 1 , ω 2 , a "moving standing wave", to induce a two-photon stimulated Raman transition through an excited state |r i (see Fig. 3). When the laser beams are far off-resonant with respect to this dipole-allowed transition, such that the detuning ∆ is much larger than the decay rate of the excited state Γ and the laser Rabi frequencies Ω 1,s , Ω 2,s , where s =↑, ↓, it is possible to eliminate the excited state from the dynamics and obtain a Hamiltonian that only involves the spins and the phonons (see Appendix A for the effective master equation). In particular, when the laser beatnote is ω L = ω 1 − ω 2 ≈ ω y ω 0 ∆ ( Fig. 2(a)), one obtains the following ion-laser Hamiltonian Figure 3. Lambda scheme for the dipole force: Two electronic states {| ↑ j , | ↓ j } with an energy difference ω 0 are selected from the ground-state manifold n 2 S 1/2 to form an effective spin-1/2. These states can be manipulated by a pair of laser beams with Rabi frequencies Ω 1,σ , Ω 2,σ that induce a transition to an excited state |r j in the n 2 P 3/2 manifold from the spin state s =↑, ↓ in n 2 S 1/2 . Note that the transition wavelength λ sp lies in the optical regime. When the beatnote ω L = ω 1 − ω 2 is tuned close to the trap frequency ω y , such that the detuning ∆ is much larger than the spontaneous decay rate Γ of the transition and the Rabi frequencies, one obtains the desired state-dependent dipole force. which can be interpreted as a time-dependent differential ac-Stark shift due to the interference of the two laser beams. Here, we have introduced σ z i = |↑ i ↑ i | − |↓ i ↓ i |, and the two-photon differential Rabi frequency Ω L = (Ω 1,↓ Ω * 2,↓ − Ω 1,↑ Ω * 2,↑ )/2∆. The corrections due to the spontaneous decay from the excited level contribute to the Rabi frequency with Ω L → Ω L (1 + (Γ/∆) 2 ), and can be thus neglected if Γ ∆.
The idea is that by controlling experimentally the polarizations, intensities and detunings of the laser beams, one can find regimes where the effective Rabi frequency Ω L does not vanish and leads to a spin-dependent dipole force. Let us note that the adiabatic elimination also leads to standard ac-Stark shifts and to a running wave that only couples to the vibrational excitations. These terms must be compensated by carefully selecting the laser-beam parameters. Let us finally emphasize that the effective Raman wavevector k L = k 1 − k 2 = k L cos θ e x + k L sin θ e y has a component along the ladder plane (see Fig. 2(b)), and in principle couples to both the planar and transverse phonons. The reason for the announced decoupling from the planar vibrational modes is that the beatnote of the laser beams, which is near the resonance of the transverse phonons, will lie far off-resonance with respect to the planar vibrational modes (Fig. 2(a)). This qualitative argument will be quantified below, and supported numerically in Sec III.
After introducing the phonon operators in Eq. (7), we perform a Taylor expansion of Eq. (10) for the small transverse and planar Lamb-Dicke parameters By setting ω L ω y , such that the bare detuning fulfills |δ y | = |ω y − ω L | |δ x | = |ω x − ω L |, ω y , one can neglect all nonresonant terms apart from a state-dependent dipole force that couples the spins to the transverse phonons. In the interaction picture with respect to the phonon Hamiltonian (9), we get where δ n⊥ = Ω ⊥ n − ω L . We note that, in order to neglect all the remaining terms of the Taylor expansion, a rotating wave approximation (RWA) must be performed, provided that The first condition is required to neglect the off-resonant contributions to the ac-Stark shift of the energy levels. Additionally, by virtue of Eq. (11), this condition ensures that the spin-phonon couplings involving higher-order powers of the transverse phonon operators can also be neglected. For the regime considered in this work, namely ω L ≈ ω y ω x ≥ ω z , it will suffice to consider Ω L /2π ≤ 0.1-1 MHz to accomplish this constraint. The second condition in (13) is necessary to avoid that the dipole force also couples the spins to the planar phonons. Hence, it characterizes the parameter regime where our previous qualitative discussion about the decoupling of the planar vibrational modes holds. The fulfillment of this condition relies on the large gap between the frequencies of the transverse and planar normal modes. Besides, by setting the laser-beam arrangement such that θ ≈ π/2, one obtains η n 1, which warrants the fulfillment of the last condition. In Sec. III, we will confirm the validity of these constraints numerically for the particular case of a zigzag ion ladder.
Depending on the spin state, the dipole force in Eq. (12) pushes the ions in opposite directions transversally to the ladder plane. Besides, the phase of this pushing force depends on the ratio of the ion equilibrium positions and the effective wavelength of the interfering laser beams. As shown below, this is precisely the parameter that shall allow us to control the anisotropy of the effective spin-spin couplings.

D. Spin models with tunable anisotropy and frustration
There are numerous situations in nature where interactions are mediated by the exchange of particles. Of particular relevance to the field of magnetism is the so-called Ruderman-Kittel-Kasuya-Yosida (RRKY) mechanism [54], whereby a Heisenberg coupling between distant nuclear spins is mediated by electrons from the conduction band of a metal. The sign of the Heisenberg couplings alternates between ferromagnetic/antiferromagnetic as a function of the ratio between the Fermi wavelength and the mutual spin distance. We have recently shown [18] that a similar phenomenon occurs for trapped-ion ladders subjected to the dipole force in Eq. (12). Instead of a Heisenberg-type interaction between the nuclear spins, one obtains a periodically-modulated Ising-type coupling between the spins formed by two electronic states of the ion. Interestingly enough, the modulation can be experimentally tailored by controlling the direction of propagation of the interfering laser beams providing the dipole force.
The Ising interaction between the effective spins of two distant ions can be understood as a consequence of the virtual phonon exchange between these ions. The dipole force (12) pushes the ions transversally, exciting thus the transverse vibrational modes. These phonon excitations, being collective, can be reabsorbed elsewhere in the ion crystal, providing a mechanism to couple the distant spins. The exact expression can be obtained by a Lang-Firsov-type transformation [15,[55][56][57][58][59] that decouples the spins from the phonons This canonical transformation leads to an effective Hamilto- where we consider a picture such that the vibrational Hamiltonian H p (9) absorbs the time-dependence of the dipole force H d (12), namely n a n ). This leads to the following representation of the total Hamiltoniañ (15) After applying the above Lang-Firsov-type transformation, and moving back to the original Schrödinger picture, we obtain the following effective Ising model The effective spin couplings have the following expression As announced previously, the phase dependence of the dipole force (12) on the ratio of the ion equilibrium positions and the effective wavelength of the light has been translated in the particular periodic modulation of the interaction strengths J eff i j ∝ cos(k L · r 0 i j ). This sign alternation is similar to that found in RKKY metals, such that the transverse phonons play the role of the conduction electrons, and the wavelength of the interfering laser beams acts as the Fermi wavelength. In the regime of interest for the ladders, κ y 1, the couplings display a dipolar decay law, where we have introduced such that η y = k L sin θ / 2mω y is the bare Lamb-Dicke parameter. From these expressions, it becomes apparent that the interactions between spins belonging to the same leg of the ion ladder (φ i j = 0) correspond to antiferromagnetic J eff i j > 0 Ising couplings. Conversely, the interactions between the spins from different legs of the ladder (φ i j = 0) can be ferromagnetic J eff i j < 0 or antiferromagnetic J eff i j > 0 depending on the laser parameters. We can thus tune the sign and magnitude of the spin-spin couplings anisotropically. Note that, even if the typical optical wavelengths are much smaller than the mutual ion distances λ L l z ≈ 1-10 µm, the angle θ can be tuned around θ ≈ π 2 so that φ i j attains any desired value φ i j ∈ [0, 2π]. Alternatively, it is also possible to maintain the laser-beam arrangement fixed, and modify the anisotropy ratio κ x = (ω z /ω x ) 2 in order to control the inter-ion distances Either of these two methods will turn out to be essential to explore the full phase diagram of the frustrated quantum spin ladders.
This anisotropy of the model leads to frustration when where P stands for the elementary plaquette of the lattice, which is a triangle in our case and (e)). Therefore, the Ising model is frustrated if there is an odd number of antiferromagnetic couplings per unit cell. We note that this standard criterion of geometric frustration [60] can be extended to situations where quantum fluctuations and frustration have the same source [61]. In our case, however, the frustration is purely classical, interpolating between the antiferromagnetic frustration, and the frustration due to competing ferromagnetic and antiferromagnetic interactions. In order to introduce quantum fluctuations, a microwave directly coupled to the atomic transition yields where we have introduced σ x i = |↑ i ↓ i | + |↓ i ↑ i |, and h plays the role of an effective transverse field of strength h due to the microwave, which is responsible for the quantum fluctuations.
Let us close this section by rewriting the effective anisotropic quantum Ising model (21) in a notation that is more appropriate for quantum spin ladders (see Fig. 4). We substitute the label of the spins i = 1, · · · , N for two new indices that account for the number of legs γ = 1, · · · , n l , and the number of spins in each of these legs i s = 1, · · · , L γ , such that ∑ γ L γ = N. Then, the Hamiltonian can be rewritten as H eff = H leg + H rung , where the leg and rung Hamiltonians are Here, we have introduced the intra-and inter-leg Ising coupling strengths, which are respectively Figure 4. Anisotropic spin-spin interactions: The legs of the trapped-ion ladder are labeled by the index γ = 1, . . . , n l , whereas the spins in each leg correspond to i s = 1, . . . , L γ . The effective spin ladder hamiltonian is composed of two terms: H leg contains single-spin terms corresponding to the transverse field h ≶ 0 (yellow-red laces) and the antiferromagnetic dipolar couplings J > 0 (yellow bonds). H rung contains the couplings between spins of different rungs, such thatJ ≶ (yellow-red bonds) depends on the mutual rung distance d γ µ . We also show the unit vectors a 1 , a 2 of the triangular lattice.
whereJ eff = J eff cosφ γ µ , such thatφ γ µ = k L cos θ d γ µ , and d γ µ is the inter-leg distance (see Fig. 4). These are the central equations of this manuscript, describing a general n l -leg quantum Ising ladder, whereby the each of the legs corresponds to an antiferromagnetic quantum Ising chain with long-range dipolar interactions. The one-dimensional chains are coupled by means of dipolar Ising pairwise interactions with a strength and sign that can be experimentally controlled. In particular, we shall be interested in a ferromagnetic coupling that competes with the intra-leg antiferromagnetic interactions, leading thus to the phenomenon of magnetic frustration.
At this point, it is worth commenting on the interesting recent proposal [62] for the simulation of any network of N interacting spins by using a linear crystal of ions. By exploiting N different laser beams individually addressed to each ion, such that the detunings and Rabi frequencies are fine tuned, it is possible to synthesize the connectivity of any desired network. Our approach is different since it exploits the specific geometry of self-organized planar ion crystals, and can be scaled to large ion ensembles straightforwardly. Additionally, it benefits from the simplicity of using a single laser-induced dipole force that lies far off-resonance from the whole vibrational branch. This contrasts the proposal in [62], where the forces lie within the vibrational branch, and resonance effects must be carefully avoided for larger ion chains where the phonon branches become denser. We note that the scheme in [62] has a higher flexibility in the simulated lattices, although the geometry of the ladders in our approach can be partially modified following the prescriptions of Sec. IV.
In the following section, we support the validity of this analytical treatment with a numerical study of the zigzag ladder.

III. A DETAILED CASE: THE ZIGZAG LADDER
A. Numerical support for the anisotropic Ising model We consider the simplest scenario where the anisotropy of the Ising model can be tested, namely, a three-ion chain in a zigzag configuration. We consider the following guiding numbers for the trap frequencies ω y /ω z = 20, ω x /ω z = 1.43, and ω z /2π = 1 MHz, although we emphasize that the scheme will equally work for different values as far as the above constraints are met. These values lead to the equilibrium positions We focus on the crucial assumption that allows us to derive the effective frustrated Hamiltonian (21), namely the possibility to neglect the pushing force on the planar modes. To address the validity of this approximation, we start from the vibrational Hamiltonian in Eq. (3), and introduce the creationannihilation operators for the local ion vibrations Even if the state-dependent dipole coupling in Eq. (10) only acts along the xy plane, the motion along the z-axis gets coupled through the Coulomb interaction (6). Therefore, we must treat the complete vibrational Hamiltonian The state-dependent dipole force can also be written in this local basis, yielding the following Hamiltonian where η α = k Lα / √ 2mω α are the bare Lamb-Dicke factors. In order to integrate numerically the Schrödinger equation for the timescales of interest t f ≈ 1/|J eff |, which are three orders of magnitude larger than the timescale set by the dynamics of the dipole force 1/ω L , it would be desirable to work in a picture that absorbs the fast time-dependence of (27). This is possible if we neglect the counter-rotating terms of Eq. (26) that correspond to phonon non-conserving processes, which is justified by a rotating-wave approximation when Then, it is possible to move to a picture where the creation-annihilation operators rotate with the laser frequency, |ψ(t) = exp(it ∑ αi ω L a † iα a iα )|ψ(t) , and the Hamiltonian that will be numerically exploredH =H 1 +H 2 becomes time-independent, namelỹ Note that this unitary transformation does not affect the spin dynamics, and is thus well-suited to study the validity of our previous derivation of the effective Ising model. In order to account for two sources of noise that are usually the experimental limiting factors for spin-oriented QSs, we include a fluctuation of the atomic resonance frequencỹ where ∆ε(t) is a stochastic process. This term may correspond to the Zeeman shift of non-shielded fluctuating magnetic fields, or to a non-compensated ac-Stark shift caused by fluctuating laser intensities. Its dynamics can be modeled as a stationary, Markovian, and Gaussian process [63], as follows where n g is a unit Gaussian random variable, and c, τ characterize the diffusion constant and the correlation time of the noise. For short correlation times τ t, one obtains an exponential damping of the coherences with a typical time T 2 = 2/cτ 2 . We set τ = 0.1T 2 , and T 2 ≈ 10 ms, which is a reasonable estimate for the observations in experiments. Note that this dephasing timescale still allows for the observation of the faster coherent spin dynamics 1/J eff i j ≈ 1 ms T 2 . We study numerically the time evolution under the total HamiltonianH = ∑ mHm in Eqs. (29)- (30), and compare it to the effective description H eff for a vanishing transverse-field in Eq. (21). We note that the following numerical simulations focus on a ground-state-cooled crystal, where the nine vibrational modes have one excitation at most. To consider the dephasing noise, we integrate over N = 10 3 different histories of the fluctuating frequency (31), and perform the statistical average of the dynamics. The effects of finite temperatures, together with other error sources, are addressed in Sec. IV. To observe a neat hallmark of the anisotropy due to the modulation of the interaction strengths J eff i j ∝ cos(k L · r 0 i j ), we study the dynamics of the spin state |ψ s = for two sets of parameters. (i) Allowed spin hopping: We set the following parameters η y = 0.1 = 10η x , ω L = 1.1ω y , and Ω L = 0.15|δ y |/η y , and direct the laser beams so that θ = π/2 (i.e. e x · k L = 0). In this case, the interfering radiation does not propagate along the triangle plane, and there is no modulation of the sign of the couplings. We observe that the initial spin excitation |+ located at site 1, can tunnel to the two remaining sites of the triangular plaquette as a consequence of the Ising-like coupling. In Fig. 5(a), we compare the numerical results with the effective description. From this figure, we can conclude that the effective Ising Hamiltonian yields an accurate description show a remarkable agreement with the numerical simulation of the complete HamiltonianH ( σ x 1 (t) red circles, σ x 2 (t) blue triangles, σ x 3 (t) yellow squares), and both describe the transport of the spin excitation around the plaquette. (b) In the inhibited hopping regime, one also observes a complete agreement between both descriptions, which describe how the spin excitation cannot occupy site 2 as a consequence of the vanishing Ising interaction. of the spin dynamics in a timescale t f ≈ 1/J eff -2/J eff ≈ 1-2 ms< T 2 ≈ 10 ms. For longer timescales, the dephasing leads to a larger deviation from the effective description. Note that this result supports the validity of the isotropic Ising interaction, but we still have to address whether our scheme to control interaction anisotropy is also accurate.
(ii) Inhibited spin hopping: We use the same parameters as before, but now consider that the interfering laser radiation also propagates along the plane defined by the triangle. In particular, we set the angle and the effective wavelength as , which leads to the factors φ 12 = φ 32 = π/2, such that the effective Ising couplings between sites 1-2 and 2-3 completely vanish J eff 12 = J eff 23 = 0. Therefore, we reach a highly anisotropic situation where the spin excitation can only hop between sites 1 ↔ 3. This is confirmed by the numerical results in Fig. 5(b), where a good agreement with the effective description is displayed once more in the timescale t f ≈ 1-2 ms. Let us finally stress that an experiment with three ions in this triangule would be a neat proof-of-principle to show that anisotropic Ising models can be realized following our scheme. We also stress that this scheme is amenable of being scaled to larger systems, since the periodically modulated cou-  plings are controlled globally by the ratio of the laser wavelength to the inter-leg equilibrium positions, and thus does not depend critically on the size of the Coulomb crystal.

B. Frustration by competing interactions
Once the validity of the effective Ising model (21) has been numerically supported, we can exploit the anisotropic interactions (17) to interpolate between a frustrated Ising ladder due to antiferromagnetic couplings, or due to the competition of ferromagnetic and antiferromagnetic interactions. In Fig. 6(a), we present a scheme of the spin frustration. The intra-leg spin couplings in the ladder (yellow) correspond to antiferromagnetic interactions which, according to Eq. (17), cannot be modulated. On the other hand, the inter-leg couplings along the diagonal rungs of the ladder may correspond to antiferromagnetic (yellow) or ferromagnetic (red) interactions depending on the value of φ i j = k L cos θ (x 0 i − x 0 j ). Both situations lead to frustration (see Fig. 6(b)), since only two of the bonds can be satisfied simultaneously. By using the normal modes of the inhomogeneous zigzag chain, we compute numerically the spin couplings, and show that for φ j 0 , j 0 +1 = 0, we obtain an antiferromagnetic coupling (Fig. 6(c)), whereas an alternating sign arises for φ j 0 , j 0 +1 = π/2 ( Fig. 6(d)).

IV. EXPERIMENTAL CONSIDERATIONS
In previous sections, we discussed the regimes of validity of the spin models by both analytic and numerical methods. Here, we analyze the capability of ion-trap experiments to meet the required conditions, considering current technology limitations and possible sources of error in ion-trap experiments. This study is supported by initial experiments.

A. Trap design study
A crucial condition for the frustrated quantum spin models is ω y ω x ≥ ω z , which leads to the clustering of vibrational branches exploited in Sec. II. This condition implies that the rotational symmetry of the trap potential in the xyplane, which is a common property of linear Paul traps, must be explicitly broken. In this section, we discuss two possible strategies to achieve this goal, and present supporting evidence based on numerical and experimental results.
The first approach can be implemented in any linear Paul trap, such as the symmetric electrode arrangement shown in the inset of Fig. 7(a). By applying a positive offset voltage U offset on the DC electrodes, the trapping potential becomes stronger along the diagonal direction joining the DC electrodes, and the trapping frequencies fulfill ω y > ω x (and vice versa for a negative offset voltage). For the experimental results shown in Fig. 7(a), resonant radio-frequency radiation was used to excite the different modes, such that the induced ion motion was observed on a CCD camera. This method allows for the estimation of the trap frequencies, which show a clear agreement with numerical ion-trajectory simulations [64](main panel of Fig. 7(a)). Note that the anisotropy of the radial frequencies ω x /ω y can be tailored by the offset voltage, which also rotates the trap axes, and effectively changes the angle θ of the laser wavevector k L with the x-axis (see Fig. 7(a)). This pinpoints the possibility of shaping the spin frustration according to Eqs. (18)- (19) by modifying the electrode voltages, avoiding thus the more demanding modification of the laser-beam arrangement. In Fig. 8, we show the measured positions for a N = 17 ion crystal. We show how the trap-frequency anisotropy can be exploited to synthesize a particular 3-leg ladder whose equilibrium positions match perfectly the numerical predictions.
The second strategy is to exploit a trap design with a nonquadrangular arrangement of the electrodes. This breaks directly the rotational symmetry of the trapping potential (see the inset of Fig. 7(b)). Both the experimental data, as measured by laser spectroscopy, and the numerics indicate a sizable splitting of the radial frequencies. Note that the dependence of the angle θ on the offset voltage depends strongly on the trap geometry. Unfortunately, in the present case, the region with the largest tunability of θ still coincides with the lowest radial anisotropies.
In order to optimize both effects, we have designed a new trap (Fig.7(c)) that meets with the special requirements of the QS. While the three-ion case-study may be realized with any state-of-the-art experimental setup, the larger scale QS with planar spin systems with about 50 ions, as shown in Fig. 1, will require a special trap design optimized according to the following conditions: (a) the trapping potential should be highly anisotropic ω y ω x , (b) the trap should allow for the confinement of a large number of ions N < 100, (c) trap frequencies should be high such that an initialization of the ion crystal in a low thermal motional state is possible, (d) the trap geometry should reduce excess micromotion, and (e) optical access should allow for readout of the spin state. We now discuss the methods to achieve these requirements.
The geometry for such a trap device is sketched in the inset of Fig. 7(c), where four round bars (r = 0.625mm) form the radial trapping potential. As we choose a non-quadrangular arrangement with distances (d x , d y ) = (0.42 mm, 2.90 mm), the two simulated radial frequencies directly become nondegenerate. With a trapping drive frequency of Ω rf /2π = 22 MHz, and amplitude of 2.5 kV, we find a sufficiently high anisotropy ω x /2π = 0.43 MHz< ω y /2π = 2.01 MHz along the requirement (a). The axial confinement is generated by two endcaps seperated by a distance of 25 mm, which lead to ω z /2π = 0.23 MHz with 2 kV applied on the endcaps. In this potential, a N = 4 ion crystal undergoes the first structural transition to the zigzag configuration, and for N > 9 the second structural transition is obtained. We have calculated the equilibrium positions and vibrational modes of a threelegged ladder of N = 19 ions for these particular trap frequencies (Fig. 9). We note that by lowering the axial DC voltage, we can increase N for different crystal structures, fulfilling the requirement (b). Let us now address condition (c). By setting the trap frequency to ω y /2π =2.01 MHz, the mean phonon numbers after Doppler cooling lie belown y < 5. In order to have a thermal error below 1% for these phonon numbers (see Eq. (C7)), the spin-phonon coupling should be smaller than Ω L = 0.02|δ y |/η y . However, this reduces the spin-spin interactions and magnetic-field noise may affect the dynamics at the corresponding long timescales. Therefore, multi-mode EIT cooling [65] ton y = 0.1 shall allow us to keep the error rate to 1%, while maintaining spin-spin interactions in the kHz-regime (i.e. Ω L = 0.15|δ y |/η y ). Let us stress that the thermal error may be minimized by considering evolution times that are multiples of the detuning of the closest vibrational mode. In order to estimate the effects of the excess micromotion according to the point (d), we have estimated the ξ i values to be ξ 2i = 0 and |ξ 1i | < 0.05π with the same beam angles as in Eq. (B10). This should lead to relative errors on the order of 4-5%. Finally, the inter-ion distances are about 10 µm for the central ions, such that we can meet with requirement (e) via high numerical aperture optics allowing for single-site readout.
With the above parameters, the relative orientation of the ion crystal and the laser wavevector k L yields an angle of θ = π/2. As discussed above, we can vary θ and tune the spin frustration by applying an offset to the DC electrodes, which rotates the crystal with respect to the fixed wavevector k L . As shown in the simulations presented in Fig. 7(c), where the new trap design allows for a smooth tunability of θ , while preserving a strong anisotropy ω x /ω y . With U offset = −4.9 V, we find ω x /ω y = 0.21, while the angle becomes θ = 0.49π,  In Figs. 1(b), (c) and (f), the ions self-organize naturally in a geometry of bond-sharing triangles. By controlling the number of legs in such triangular ladders, the QS is already capable of exploring a variety of interesting cooperative phenomena (see Secs. V and VI). However, it would be highly desirable to have a method to modify these geometries, widening thus the applicability of our quantum simulator.
Since the geometry of the spin model (21) is determined by the couplings J eff i j between the spins σ z i ↔ σ z j , a possibility to tailor the geometry is to switch on/off some of these interactions. A possible route that allows for such a selective coupling is a type of laser-beam hideout. Thanks to the relatively large distances between the ions l z ≈1-10 µm, it is possible to address them individually with laser light [98] (note that the tools used so far only require global addressing). The idea is that, after the global optical pumping to the state |ψ s = ⊗ i | ↓ i , focused laser radiation will selectively transfer the population to a different level |h i that is not coupled to the spin-dependent force. Therefore, the ion gets effectively hidden. This may be achieved by polarization selection rules, or alternatively, by the ac-Stark shift of highly-detuned laser beams. Using this idea, it is possible to construct a simple ladder of corner-sharing triangles (see Fig. 1(c)), or a stripe of the two-dimensional Kagome lattice (see Fig. 1(d)). Additionally, it opens the possibility of introducing defects in the lattice in order to study the role of disorder. As discussed in Sec. V, this tool opens many possibilities for our QS. C. Imperfections and noise in the quantum simulator Let us now address the possible sources of error in the QS. We place a special emphasis on the ion micromotion, which has not been discussed in previous QSs [8][9][10], since it can be cancelled for linear chains. We also discuss other error sources shared with the linear-chain QS, such as thermal fluctuations and heating of the phonons, dephasing of the spins, and photon scattering. We stress that, provided that the following constraints to minimize the micromotion errors are fulfilled, the proposed ladder QS should not present more limitations than its linear-chain counterparts [8][9][10].
(i) Sources of dephasing: Usually, the most important important source of error in the experiments is the dephasing of the electronic states. These might be caused by fluctuating Zeeman shifts on the magnetic-field sensitive states required to implement the spin-dependent dipole force (10); or by fluctuations in the laser-beam intensities leading to uncompensated ac-Stark shifts. A reasonable estimate of the typical coherence times T 2 ≈ 5-10 ms [6] shows that these terms are important error sources in the timescales of interest t f ∝ 1/J eff , where J eff /2π ≈ 1 kHz. As shown in Fig. 5, the desired coherent dynamics described by the effective Ising model (21) dominates the behavior of the system for t f ≈ 1-2 ms. This sets the timescale for the QS of frustrated magnetism to be in the millisecond range, after which the dephasing is so strong that the quantum simulator is no longer faithful.
An interesting protocol to overcome both sources of dephasing simultaneously is based on the concept of continuous dynamical decoupling already applied to trapped ions [73]. In our case, it could be implemented by a strong driving of the carrier transition [74]. Hence, instead of relying on the spin-dependent force (10), it suffices to combine a red-sideband term with the strong driving of the carrier. In order to obtain a quantum Ising model, one should introduce an additional microwave that provides a Zeeman shift oscillating at the Rabi frequency of the strong carrier driving. The phase diagram can be explored adiabatically by an intermediate spin-echo pulse that refocuses the fast oscillations due to the strong driving.
(ii) Role of the ion micromotion: The initial Hamiltonian in Eq. (1), upon which the derivation of the effective quantum spin ladders (22) has been built, is based on the so-called pseudo-potential approximation [6]. This approximation neglects the effects of the ion micromotion: a fast motion of frequency Ω rf that is synchronous to the radio-frequency (r.f.) field of the Paul trap. This allowed us to focus directly on the slow secular motion of the trapped ions, which is described by an effective harmonic potential (1) with frequencies ω α Ω rf . The validity of this assumption is well justified for single trapped ions, where there are experimental methods to minimize the effects of the micromotion [6]. For linear ion chains, the micromotion is also minimized when the equilibrium positions lie along the trap axis. Conversely, for the ion ladders considered in this work, the equilibrium positions necessarily lie off the trapping axis, leading to an additional micromotion that cannot be compensated.
In Appendix B, we present a detailed discussion of the conditions under which this micromotion does not affect the QS. Let us comment on the possible sources of imperfections, and the conditions to minimize them. The r.f. heating of the transverse phonon modes responsible of the spin couplings is negligible when a condition easily verified for the parameters in this work. In addition, during the cooling stage of the transverse vibrational modes, the laser frequency must be carefully tuned to avoid heating caused by the micromotion-induced broadening of the transition, or to additional micromotion sidebands [72]. We have also considered how the micromotion may affect the spin-dependent dipole force (10) and the spin-phonon coupling (12). To neglect the associated contributions, it is necessary to consider that the laser beams lie almost parallel to the y-axis, and are far off-resonant (see Fig. 3), namely ω L /2π ≈ ω y /2π ≈10 MHz Ω rf /2π ≈ 0.1 GHz ∆/2π ≈ 10 GHz. Besides, the Rabi frequencies of the laser beams should fulfill Finally, some care must also be placed in order not to induce two-photon transitions to different states |a i of the atomic ground-state manifold. To avoid these processes, the relative Zeeman shifts of such transitions δ a,s , and the associated Rabi frequencies must be controlled such that |Ω 1,a Ω * 2,s | |δ a,s − ω L ± Ω rf |.
We thus conclude that, provided that the above restrictions on the parameters are fulfilled, the unavoidable excess micromotion of the ions in a ladder geometry should not affect our QS.
(iii) Thermal motion and heating of the ions: In Sec. III, we have supported the validity of the QS based on numerical results for ground-state cooled ion crystals. However, ever since the early schemes for phonon-mediated gates, the thermal fluctuations of the ions have been identified as a potential source of errors that must be carefully considered.
In Appendix C, we estimate the thermal error of our QS by both analytic and numerical methods. We solve exactly the Heisenberg equations of motion for the spin-dependent dipole force (15), and derive a scaling law for the thermal contribution to the T = 0 spin dynamics, which is then supported numerically. We find that the relative thermal error, defined as wheren y represents the mean number of transverse phonons in the center-of-mass mode. By considering the parameters used in section III, the thermal fluctuations have a small contribution (1%) forn y ≈ 0.1 (see Fig. 19), which shows that perfect ground-state cooling is not required to implement the QS accurately. At this point, we should mention that the above arguments are valid for a vanishing transverse field h = 0. We note, however, that the contribution to the thermal error due to h > 0 in the regime of interest h ≈ J eff has been shown [15] to be negligible with respect to the estimate in Eq. (35).
In Appendix C, we also present a phenomenological master equation that models possible heating mechanisms in the ion trap. We show that the relative heating error, defined as where Γ h is the heating rate, only yields a small contribution (1%) to the overall error of the QS for heating times above 5 ms/ phonon (see Fig. 20).
(iv) Spontaneous photon scattering: The spontaneous emission can severely limit the advantage of quantum-information protocols [75][76][77]. In our case, the use of two electronic ground-states as the effective spins (see Fig. 3) makes the direct spontaneous emission | ↑ i | ↓ i negligible. However, due to the Lambda-scheme responsible for the dipole force, there can be photon scattering events from the intermediate excited state |r i , which must be carefully considered.
In Appendix A, we describe in detail the derivation of an effective master equation where H d corresponds to the dipole force (10), and accounts for the two possible decoherence channels. These are the so-called Raman and Rayleigh scattering of photons, and are contained in the following effective jump operators L eff 1 = √ Γ ∆ Ω 1,↓ e i(k 1 ·r−ω 1 t) + Ω 2,↓ e i(k 2 ·r−ω 2 t) |↓ ↓|+ + √ Γ ∆ Ω 1,↑ e i(k 1 ·r−ω 1 t) + Ω 2,↑ e i(k 2 ·r−ω 2 t) |↓ ↑|, where we have assumed that the scattering rate is much smaller than the laser detuning Γ ∆. Note that the first term of each jump operator does not change the spin state (i.e. Rayleigh scattering), being responsible for pure dephasing. Conversely, the second term flips the spin state (i.e. Raman scattering), leading thus to damping of the spin populations.
A conservative estimate of the photon scattering is to consider the individual effective rates (37), which scale as Γ eff = Γ(|Ω l,s |/∆) 2 ≈ |Ω L |(Γ/∆). By using the parameters introduced in Sec. III, we find that Γ eff /J eff ≈ 10 3 (Γ/∆). Hence, by considering a sufficiently large detuning, it is possible to keep the scattering rates below Γ eff /J eff < 10 −1 , so that they do not compromise the accuracy of the effective description (21). For instance, for typical decay rates Γ/2π ≈1-10 MHz, one must consider detunings in the range ∆/2π ≈10-100 GHz.
(v) Spatial dependence of the laser-beam profile: The effective spin couplings in Eq. (19) assume a Rabi frequency that is constant along the ion crystal. For large crystals, however, the characteristic gaussian profile for the laser intensities may lead to weaker Rabi frequencies on the boundaries of the crystal. This effect will give rise to inhomogeneous spin-spin couplings that must be added to the inhomogeneities caused by the varying inter-ion distance in Coulomb crystals. Rather than considering these terms as an error, they can be seen as a gadget that makes the many-body problem even more interesting. In particular, they will lead to inhomogeneous critical points which may be responsible of interesting effects [46].

D. Efficient detection methods
A crucial part of a QS is the ability to perform measurements that yield information about the Hamiltonian under study. For the frustrated quantum Ising ladders (22), a QS would start by preparing the so-called paramagnetic state |P = ⊗ i |→ i with |→ i = (|↑ i + |↓ i )/ √ 2, which is the ground-state of the model in the absence of spin interactions. Such a separable state can be accurately prepared by means of optical pumping, followed by π 2 -pulses globally addressed to the whole ion crystal [6]. This step is followed by the adiabatic modification of the Hamiltonian parameters J eff ,J eff , h, in order to connect the paramagnet to other phases of the quantum Ising ladder (see Secs. V and VI below).
A direct approach to the measurement of these phases is the so-called quantum state tomography, more precisely, the full determination of the state of the system [67]. However, this approach becomes highly inefficient for many-body systems due to the exponential growth of the composite Hilbert space, and alternative schemes must be studied. An interesting possibility for state estimation are the methods based on matrixproduct representations of the states [68]. Another alternative that does not require full quantum state tomography is the measurement of order parameters characterizing the phases.
One of the advantages of trapped-ion experiments with respect to other platforms is their ability to perform highlyaccurate measurements at the single-particle level [6]. The technique of state-dependent fluorescence allows for the measurement of single and joint probability distributions of the electronic states P ↑ i , P ↑↑ i j . From the spatially resolved fluorescence, it is possible to infer local expectation values, such as the magnetization σ z i = 1 2 (P ↑ i − 1), or two-body correlators As discussed in Sec. VI, these observables usually contain all the relevant information about the different phases. In particular, one could study the dependence of the correlator with the distance (48), or infer the magnetic structure factor S zz (q) = ∑ i j σ z i σ z j e iq(i− j) . Let us now comment on the possibility of recovering some of these magnitudes from global properties of the fluorescence spectrum. By measuring the probability to find a fraction of n-ions in the excited state P ↑ (n) [10], one obtains the total magnetization without single-site resolution m z = 1 N ∑ i σ z i . To measure correlators, note that the fluorescence of an ensemble of emitters may carry information about their correlations [69]. In our case, the resonance fluorescence associated to the cycling transition depends on the collective properties of the ion ensemble [70]. In fact, the power spectrum in a particular detection direction,r, is related to the structure factor (38) where λ is the wavelength of the emitted light. Even if the ion spacing is much larger than the optical wavelength l z λ , one may compensate it by setting the photodetector almost orthogonal to the plane defined by the ladder, and by using photodetectors with a very good angular resolution. We finally note that such magnetic structure factors yield a lower bound on the entanglement without the need of state tomography [71].

V. SCOPE OF THE QUANTUM SIMULATOR
Once the validity of the spin-ladder Hamiltonian (22) has been addressed by both analytic and numerical methods (Secs. II and III), and its experimental viability discussed (Sec. IV), we can now focus on the many-body models to be explored with the QS. We place a special emphasis on the range of collective phenomena that are not fully understood, or have not been addressed so far to the best of our knowledge. These would directly benefit from the advent of such a QS.
A. J 1 -J 2 quantum Ising model The first many-body model that may be targeted with the proposed QS is the so-called axial next-to-nearest neighbor Ising model [78], supplemented by quantum fluctuations (see [34] and references therein). It has the Hamiltonian which consists of nearest neighbor (J 1 < 0) and next-tonearest neighbor (J 2 > 0) couplings, and the transverse field h.
In spite of the mild looking appearance of this Hamiltonian, it has a rich phase diagram with some features that are still controversial. We refer to this model as the J 1 -J 2 quantum Ising model (J 1 -J 2 QIM), a paradigm of frustrated magnetism. The Hamiltonian of our QS (22) corresponds to the J 1 -J 2 QIM for the simplest possible geometry, namely, the twoleg zigzag ladder [18]. The indices of this ladder, γ = 1, 2, i s = 1, · · · , N/2 (Fig. 4), are mapped onto a one-dimensional chain by the relation i = 2(i s − 1) + γ (Fig. 5(a)). Then, the analogy with the J 1 -J 2 QIM follows directly provided that Note that due to the particular laser-beam arrangement presented in Sec. II, it is possible to tailor the ratios J 2 /J 1 and h/J 1 experimentally. We emphasize that, even if a solid-state material is found to be described by such a model, a similar microscopic control of the couplings seems rather difficult to achieve. Therefore, the trapped-ion platform is ideal to explore the different regions of the phase diagram [79]. Of particular relevance is the region around the frustration point f c = J 2 /|J 1 | = 1/2, which is characterized by a macroscopicallydegenerate ground-state. Therefore, even a small amount of quantum fluctuations due to the transverse field may lift the classical degeneracy leading to a variety of magnetic phases. Such a rich phase diagram is studied numerically in Sec. VI, where we identify some additional features caused by the dipolar range of interactions present in the trapped-ion QS. Let us briefly mention that these long-range interactions introduce incompatible sources of frustration capable of stabilizing a new order that complements the ferromagnetic, dimerized antiferromagnetic, and floating phases that are also present in the short-range model (39). Hence, our QS will be of the utmost interest to explore the interplay between frustration, quantum fluctuations, and long-range interactions. Besides, the proposed QS shall be able to address some open questions about the phase diagram that are still a subject of controversy [34], such as the extent of the floating phase and the existence of a multi-critical Lifshitz point.

B. Dimensional crossover and quantum dimer models
An ambitious enterprise is the understanding of the dimensional crossover from the two-leg quantum Ising ladder onto the two-dimensional (2D) triangular quantum Ising model where the spins are labelled according to the 2D Bravais lattice vectors r m,n = ma 1 + na 2 (see Fig. 4), and we consider anisotropic couplings J 1 < 0 and J 2 > 0. The correspondence with our trapped-ion QS (22) is straightforward if one considers that m labels the spins within each leg of the ladder coupled by J 2 ↔ J, and n the different legs coupled by J 1 ↔J. Note that, following recent experimental efforts [12,80], the triangular QIM may also be realized with ions in Penning traps [81]. However, the study of the ladders and the crossover phenomena seems to be better suited to ions in Paul traps.
The triangular classical Ising model can be considered as the backbone of frustrated magnetism [32]. Already in the absence of quantum fluctuations, there are suggestive questions that deserve a careful consideration. For instance, the frustration point of the two-leg zigzag ladder f c = 1/2 flows to the isotropic point f c = 1 in the 2D model. Moreover, the macroscopic ground-state degeneracies of these models yield different ground-state entropies. We believe that it would be fascinating to explore these topics with trapped ions, which allow for the consecutive increase of the number legs (Figs. 1(c),(f)).
The dimensional crossover is even more exotic when quantum fluctuations are included. Let us note that the dimerized antiferromagnet of the two-leg zigzag ladder corresponds to two possible classical dimer coverings of the ladder (see Sec. VI), where each dimer corresponds to a nearest-neighbor bond that is not satisfied due to the frustration. Since the ground-state tries to minimize the number of dimers, each site of the dual lattice belongs to only one dimer, and the covering consists of the dimer arrangement along the rungs of the ladder. Note that each spin is connected to three satisfied bonds, and only one broken bond. The situation gets more interesting for the 2D quantum Ising model, since there, a single spin may be connected to the same number of satisfied and broken bonds. In this case, the transverse field can flip the spin and produce a resonating effect for neighboring dimers [37], providing a beautiful connection to the so-called quantum dimer models (see [82] and references therein).
Quantum dimer models were introduced [83] in the context oh high-temperature superconductivity. Here, these models provided a neat playground where to study the quantum spin liquid phases of the undoped cuprates, which were conjectured to play a key role in the onset of superconductivity upon doping [24]. However, they have evolved into an independent subject displaying exotic effects, such as topologically ordered phases and fractional excitations. Although originally introduced for Heisenberg antiferromagnets, where the dimers correspond to spin singlets, there is also a link to frustrated quantum Ising models by mapping the dimers to the broken magnetic bonds due to the frustration [37].
Their connection to the isotropic triangular quantum Ising model [37] has allowed to predict a quantum version of the phenomenon of 'order by disorder' [84], which yields an ordered phase out of the classically disordered frustration point f c = 1 when quantum fluctuations are switched on. In contrast, by switching the transverse field on the frustration point f c = 1/2 of the two-leg zigzag ladder, one only obtains a disorder paramagnetic phase (see Sec. VI). Therefore, the trapped-ion QS offers a unique playground to study how the phenomenon of 'order by disorder' sets in as the number of legs is increased, and how the anisotropy and the long range of the interactions affects it.

C. Quantum spin liquid phases
In the models considered so far, strong quantum fluctuations trigger a phase transition connecting an ordered phase to the uninteresting disordered paramagnet. However, this behavior does not exhaust all possibilities. Quantum fluctuations may be responsible of stabilizing more exotic phases that do not break any symmetry of the Hamiltonian, the so-called quantum spin liquids [85].
We have discussed how the two-leg zigzag ladder consisting of bond-sharing triangles allows for the QS of the paradigm of FQIM. However, this model only accounts for a disordered paramagnet. There exists another simple two-leg ladder which, although not so widely known, has been argued to provide the simplest instance of a frustration-induced quantum spin liquid [86] (see also [37]). This is the so-called sawtooth quantum Ising model, which consists of corner-sharing triangles described by the following Hamiltonian Here, J 1 < 0 represents the coupling along the diagonal rungs, and J 2 > 0 stands for the interactions along the lower leg of the ladder. In the isotropic frustration point f c = 1, the ground state is found by minimizing the number of broken bonds in each triangle independently, which leads to a macroscopic ground-state degeneracy. The numerical results in [86], which are based on the exact diagonalization of small ladders and perturbative expansions, support the absence of any symmetry breaking as the transverse field is increased. This effect has been coined as 'disorder by disorder', and would provide a testbed for a disordered quantum spin liquid state [37]. In order to perform a quantum simulation of this model, the geometry of a three-leg triangular ladder must be modified according to the method presented in Sec. IV (see Fig. 1(d)). Then, the model Hamiltonian (42) would follow directly from the QS Hamiltonian (22) with the usual identifications J 2 ↔ J, and J 1 ↔J. In addition to the possibility of reaching larger spin ladders to test the conclusions of [86], the trapped-ion QS may explore of the effects of f 2 = J 2 /|J 1 | = f c and longrange interactions, hopefully leading to a rich phase diagram which, to the best of our knowledge, still remains unexplored.
At this point, we should remark that quantum spin liquid phases can also be realized in (quasi) one-dimensional Heisenberg magnets [22,26]. A much harder task is to find their higher-dimensional counterparts [85]. By using the same technique to modify the geometry of the trapped-ion ladder, one may construct three-leg ladders such as the one displayed in Fig. 1(e). Note that this amounts to a single stripe of the well-known two-dimensional (2D) Kagome lattice. Hence, our QS allows for the exploration of the dimensional crossover towards the Kagome quantum Ising model. In this case, the interplay between frustration and quantum fluctuations has been argued to give rise to a 2D quantum spin liquid phase [37], which could be targeted with the proposed QS. Besides, the study of the dimensional crossover and the long-range interactions is likely to introduce a variety of interesting effects.
Before closing this section, let us remark that in addition to the aforementioned static phenomena, the trapped-ion QS is also capable of addressing dynamical many-body effects. In particular, the microscopic parameters of the above Hamiltonians can be tuned dynamically across the different quantum phase transitions. The breaking of the adiabatic approximation associated to a quantum critical point has been the subject of recent interest for different quantum systems (see e.g. [87]).

VI. A DETAILED CASE: THE J 1 -J 2 QUANTUM ISING MODEL
In this Section, we focus on the many-body physics of the J 1 -J 2 QIM. For the sake of completeness, we first review the properties of the short-range model, and then discuss the changes introduced by the dipolar range of the interactions, as realized in the trapped-ion quantum simulator.
According to the Toulouse-Villain criterion (20), this leads to a frustrated quantum spin model F = −1, whereby the ratio f controls the frustration, and g the quantum fluctuations.
The classical J 1 -J 2 Ising chain, obtained by setting g = 0 in the above Hamiltonian, can be solved exactly (see [88] and references therein). Such a solution yields two possible phases. For f 2 < 1 2 , one lies in a ferromagnetic (F) phase whose ground-state manifold is two-fold degenerate where the spins have been arranged according to the trappedion zigzag layout. Conversely, for f 2 > 1 2 , the ground-state is a dimerized antiferromagnet (dAF) with four-fold degeneracy Note that the degeneracy of the ferromagnetic manifold is related to the global spin-flipping Z 2 symmetry of the Hamiltonian U = i σ x i , such that Z 2 = {I,U} is the smallest cyclic Abelian group. On the other hand, the doubling of the dimerized-antiferromagnet degeneracy is accidental (i.e. not related to symmetries). Such an accidental degeneracy becomes more important at the point f c = 1 2 , where the frustration leads to a macroscopically degenerate ground-state. In fact, the degeneracy has been shown to scale exponentially with the number of spins d c ∝ ϕ N , where ϕ = 1 2 (1 + √ 5) is the golden ratio [89]. This yields the hallmark of frustrated magnets, namely, a non-vanishing ground-state entropy.
We are interested on the impact that quantum fluctuations may have on these degeneracies. For g 1 f 2 , the ground state corresponds to a single paramagnetic (P) state with all spins pointing towards the direction of the transverse field where |→ = (|↑ +|↓ )/ √ 2. The situation gets more interesting for intermediate fields, whereby additional exotic phases and a variety of quantum phase transitions occur. Since the model is no longer integrable, the analysis of the full phase diagram has been a big challenge, requiring the combination of a variety of techniques. For instance, the mapping to a classical 2D model [90] yields a direct link to commensurateincommensurate thermal phase transitions [88,91]. Hence, it is possible to use the methods developed in this area in order to understand the magnetic phases of the quantum model, being numerical Monte Carlo [92] and free-fermion approximations [93] two representative examples. Together with the more recent application of bosonization techniques [94] and numerical renormalization group methods [95], these studies yield the rich phase diagram represented in Fig. 10. Note that in addition to the aforementioned phases, there is a modulated paramagnetic phase (mP), and a highly-debated incommensurate floating phase (FP). Since all these different phases meet at f c = 1 2 , this macroscopically-degenerate point is also known as the multi-phase point. There are several critical lines emerging from the multi-phase point, which give rise to quantum phase transitions of second order, Kosterlitz-Thouless [96], or Pokrovski-Talapov [97] type.
Notwithstanding these big efforts, we emphasize that there still exists some controversy about the floating phase. In particular, the extent of the floating phase is still a question of debate. Whereas some results point towards a FP that prolongs towards f 2 1, other treatments predict a finite region that terminates in a multi-critical Lifshitz point (see [34] and references therein). This makes a QS of the utmost interest to settle down these discrepancies. Moreover, as we discuss in detail below, the introduction of long-range interactions leads to additional open questions that have not been previously addressed to the best of our knowledge. From our numerical survey, we conjecture that the dipolar-range of the couplings leads to the splitting of the multi-phase point and the appearance of an intermediate phase (see Fig. 11).
B. Dipolar-range J 1 -J 2 quantum Ising model Trapped ions are an ideal platform to test the effects of longrange Ising interactions. For non-frustrated systems, even if the model belongs to the same universality class as the nearest-neighbor case, these long-range interactions may shift the critical point, and favor long-distance quantum correlations [15]. For the frustrated systems under study, the effect of long-range interactions is expected to be more significant.
By considering the dipolar range of the trapped-ion ladder Hamiltonian (22), the J 1 -J 2 QIM (43) must be modified to where we have introduced the ratios f δ = J δ /|J 1 |, such that J δ account for the dipolar decay of the interactions. In this notation, the frustration ratio of the nearest-neighbor chain is f 2 = J 2 /|J 1 |, and we can study the additional frustration coming from the dipolar range f 3 , f 4 .... We have analyzed numerically the phase diagram of such a long-range frustrated spin model by means of an optimized Lanczos algorithm that allows us to reach efficiently ladders with L = 24 spins. All the results showed below have been performed including the f 3 and f 4 terms of the dipolar tail, in addition to the competing nearest-neighbor (J 1 ) and next-to-nearest-neighbor (J 2 ) couplings. We have checked that the effect of including longer ranged interactions (i.e., those corresponding to f 5 and f 6 couplings) does not affect qualitatively the results. In order to distinguish the phases, we focus on the two-body correlators, which show the following scaling when |i− j| 1 where 0 < m 0 < 1, ξ P , ξ mP stand for the correlation lengths of the paramagnets, η > 0 characterizes the algebraic decay of the correlations in the gapless floating phase, and the different modulation parameters q [... ] have been listed in Table I.   Table I. Magnetic modulation parameters q F q dAF q P q mP q FP q ? 0 According to these expressions, the F and d-AF display long-range magnetic order with different periodicities, whereas the P and mP phases are disordered. Finally, the FP phase has quasi-long range order with a modulation parameter that flows with the ratios f 2 , g (hence the adjective floating), and is generally incommensurate with the underlying lattice.
An observable capable of capturing the periodic modulations of the long-range ordered phases, and thus the features of the phase diagram, is the so-called magnetic structure factor. It is defined as the Fourier transform of the spin correlations with q ∈ [0, 2π), and should attain a maximum at the different values shown in Table I for each of the phases. In order to evaluate the magnetic structure factor numerically, we have considered periodic boundary conditions, which shall capture the bulk properties in the center of the trapped-ion ladders.
According to the scheme depicted in fig. 11, the main features of phase diagram of the dipolar J 1 -J 2 QIM can be divided into three different regions: the region where the ferromagnetic order prevails ( f 2 < 0.5), the region with the dimerized-AF order ( f 2 > 0.6), and the new intermediate region (roughly, f 2 ∈ (0.5, 0.6)) that appears due to the competition of different long-ranged frustration mechanisms. In Figs. 12(a)-(b), we show how the dAF phase is destabilized by these competing mechanisms, which accounts schematically for the splitting of the multi-phase point. We will show below how the numerical evaluation of the structure factor supports this division. We note that the structure factor has been normalized to unity.
In Fig. 13(a), we plot the order parameter S zz ( π 2 ) corresponding to the dAF phase. It can be observed that both phase transitions, namely the discontinuous transition (Pokrovski-Talapov type) separating the dAF phase from the FP, and the continuous transition (Kosterlitz-Thouless type) between the FP and the mP phases, reflect themselves clearly in the structure factor as two consecutive jumps with increasing transverse field. The fact that only one transition is observed for lower values of f 2 hint at the possibility that the multicritical point has been shifted from g = 0 to some finite value. Whether this shift survives for larger lattices, or corresponds to a finite-size effect, is an open question that cannot be addressed due to the limitations of the diagonalization routine.
The dependence of the structure factor with the lattice momentum q can be seen in Fig. 13(b). Here, we represent the evolution of this magnitude as we increase g for a fixed value of f 2 . We have chosen some representatives for each phase: for low fields g, we are well into the dAF phase, whose periodicity is given by a single q = π 2 component. Increasing the transverse field g, the state undergoes a transition to the FP, showing a ( f 2 -, g-dependent) incommensurate modulation. Upon further increase of the field, we reach the mP phase. Note that for moderate g, the modulation of this phase is apparent in the peaked form of its structure factor. Increasing further the magnetic field polarizes the spins in the transversal direction yielding a flat vanishing structure factor.
The same analysis has been carried out in the ferromagnetic region. In Fig. 14(a), we have plotted the S zz (0) component characteristic of a ferromagnet. In this case, however, only one clear phase transition arises as we increase the transverse field g for a fixed value of f 2 . Also in accordance with this observation, the modulation of the ground-state in Fig. 14(b), shows a shift between a pure ferromagnetic order for low values of g, to a modulated paramagnetic one with increasing field, and finally to the completely polarized paramagnet for g 1. This suggests the possibility that the unmodulated paramagnetic phase observed in the J 1 -J 2 QIM is not stable upon the effect of further dipolar couplings, or that its extent in the phase diagram is small enough to prevent to be accurately captured with finite-sized lattices.
In the J 1 -J 2 QIM, there exists a one-point boundary between the ferromagnetic and dimerized-antiferromagnetic phases precisely located at the muticritical point (g = 0, f 2 = 0.5). Interestingly enough, in the dipolar J 1 -J 2 QIM this is no longer true and both phases are separated by a finite region. In Fig. 15(a), we have plotted the structure factor S zz (q) for the lattice size L = 16 along a fixed value of f 2 in between the F-and dAF-phases. In this graph, a new type of modulation shows up in the form of two diferentiated peaks with momenta q 1 = π 4 and q 2 = 3π 4 . However, the precision in the determination of such peaks is limited by the size of our lattice. We have carried out additional numerics for L = 24, which allow us to bound these peaks between q 1 ∈ [ π 6 , π 3 ] and q 2 ∈ [ 2π 3 , π]. We have checked that the relative amplitude of these peaks does not depend on the ratio between the intra-and inter-leg couplings. In Fig. 15(b), we plot the amplitude of these two characteristic peaks with increasing g. The shaded area represents the range where both peaks coexist, i.e, the extent of the new ordered state. It would be very interesting to explore the precise origin of these modulations via the trapped-ion QS.
We remark that the modulations q 1 and q 2 coexist for low transverse fields. As we increase g, only the modulation with lower lattice momentum survives as a differentiated peak. Indeed, the momentum of this surviving modulation is compatible with that found in the mP phase for intermediate fields (compare the light-blue curve in Fig. 15(a) with the red ones in Figs. 14(b) and 13(b)). Whether this is a crossover or a quantum phase transition between the new conjectured phase and the mP, goes beyond the limitations of our numerical tools. All the above results have been computed in ladders with lattice parameters d = a . Note that a ladder composed of equilateral triangles would correspond to d = √ 3a. The extent of the new intermediate phase in the dipolar J 1 -J 2 QIM is indeed strongly dependent on the geometry of the triangular plaquettes. A straightforward way of measuring the extent of this phase is computing the distance between the F and dAF phases along the line g = 0. In Fig. 16, we show the order parameters S zz (0), S zz (π/4) and S zz (π/2) for different values of the anisotropy ratio d/a, which correspond to the F phase, the new intermediate phase, and the dAF phase, respectively. From these graphs (qualitatively similar results are obtained with L = 24), it is apparent that the splitting of the multiphase point is strongly enhanced in anisotropic lattices.

VII. CONCLUSIONS AND OUTLOOK
Low-dimensional quantum Ising magnets are often representative models for a variety of emergent cooperative phenomena. Unfortunately, in contrast to their Heisenberg counterparts, the identification of materials accurately described by these models has turned out to be a much more difficult task.  Moreover, whenever the models include anisotropic interactions and tunable ladder structures, the prospects of realizing them diminish even further. A promising alternative to overcome such difficulties are the so-called quantum simulators.
In this work, we have shown that cold ion crystals are promising candidates for the quantum simulation of a variety of quantum spin-ladder compounds. This avenue of research will allow for the study of collective phenomena due to the combination of frustration, quantum fluctuations, long-range interactions, and dimensional crossover. In particular, these ion crystals can be used to explore paradigmatic, yet controversial, models of quantum Ising ladders whereby analytical and numerical techniques seem to disagree. Moreover, they may also allow for the study of previously unexplored features of the models, such as the effects of long-range interactions.
First, we have shown how to control the geometry of different self-assembled trapped-ion ladders. Based on this possibility, we have presented a protocol to tailor the anisotropy of the magnetic interaction mediated by the transverse phonons, which couple pairs of distant ions indirectly. When these ions belong to different legs of the ladder, it is possible to tune both the sign and the magnitude of the spin-spin couplings by manipulating the laser-beam arrangement. The validity of this technique is supported by a detailed discussion of the possible sources of error in current ion-trap experiments, and by numerics showing an excellent agreement with our predictions.
This tool opens a vast amount of possibilities for trappedion-based QS of cooperative magnetic phenomena. For instance, we have presented a thorough description of the QS for the cornerstone of frustrated quantum Ising magnets, the J 1 -J 2 quantum Ising model. Moreover, we have also discussed how the QS has the potential of realizing quantum Ising ladders with connections to the exotic quantum dimer models introduced in the context of high-temperature superconductivity, and address the intriguing dimensional crossover phenomena. Finally, we have also pointed out how this QS may yield a route towards the long-sought quantum spin liquid phases. In this Appendix, we present a detailed derivation of the effective laser-ion interaction [50,51] in Eq. (10). By taking into account the spontaneous decay from the excited level (see Fig. 3), we can discuss the regime where the spin-dependent dipole forces arise, and also analyze the sources of error due to photon scattering in the experiments. We consider the master equation for the Lambda scheme in Fig. 3, namely The coherent part of the evolution is given by 2 Ω l,s |r s|e −iω l t + H.c., (A2) where we have introduced the energies of the internal states ε r , ε ↑ , ε ↓ , the Rabi frequencies of the transitions Ω l,s , and the laser frequencies ω l . The dissipator describing the spontaneous decay from the excited state is of the Lindblad form with the following jump operators L 1 = √ Γ|↓ r|, and L 2 = √ Γ|↑ r|, where Γ is the spontaneous decay rate from the excited state back to the spin manifold (see Fig. 3). Let us now define the detunings for all possible transitions As announced in the main text, when these detunings are much larger than the Rabi frequencies and the decay rate, namely δ l,s Ω l,s , Γ, it is possible to adiabatically eliminate the excited state from the dynamics, and obtain an effective master equation within the spin manifold. We use the formalism introduced in [99], and obtain the master equation where the Hamiltonian is H eff = H s + H r + H d , such that includes the following ac-Stark shifts and we have introduced the sum of the two decay rates Γ t = 2Γ. Additionally, we get the following two-photon stimulated Raman transitions where we have introduced σ − = |↓ ↑| = (σ + ) † , and the effective Rabi frequencies for the different Raman transitions which involve the absorption of a photon from the beam l = 1, 2 and posterior photon emission into the beam l = 1, 2. The last term of the effective Hamiltonian is the one responsible for the spin-dependent dipole forces where we have introduced the following Rabi frequencies Now, depending on the particular laser frequencies ω l , it is possible to select whether the laser-ion interaction leads to a stimulated Raman transition (i.e. ω L := ω 1 − ω 2 ≈ ε ↑ − ε ↓ =: ω 0 ), or to a spin-dependent dipole force (i.e. ω L ω 0 ). The effects of the spin-dependent dipole force are more transparent by rewriting Eq. (A10) as follows where we have introduced σ z = |↑ ↑| − |↓ ↓|, and Once the coherent part of the effective master equation has been derived, one must obtain the effective jump operators. In our setup, they can be expressed as follows Note that the first term of each of the effective jump operators corresponds to the so-called Rayleigh photon scattering, which takes place without modifying the internal spin state. In this formulation, it becomes clear why the Rayleigh scattering will only introduce dephasing when the amplitudes (i.e. terms between brackets) of each jump operator are different, as observed in recent experiments [77]. The second term of each jump operator corresponds to the Raman scattering, whereby the spin state is changed after the emission of the photon. In order to study the accuracy of this effective description (A5), we confront it with the exact numerical integration of the original master equation (A1). In Fig. 17(a), we compare both predictions for ε ↑ /ε r = 0.1, ε ↓ /ε r = 0.05, and setting the laser parameters such that δ 1,↑ /ε r = 0.5, ω L = ω 0 , Ω 2,↑ /ε r = Ω 1,↓ /ε r = 0.05, and Ω 1,↑ = Ω 2,↓ = 0. This set of parameters leads to the regime of stimulated two-photon transitions, so that we expect to find periodic Rabi oscillations in the populations with a frequency Ω eff /ε r = 2.5 · 10 −3 when the initial state is |ψ 0 = |↓ . As shown in the upper panel of Fig. 17(a), these oscillations get damped due to the spontaneous decay Γ/ε r = 0.05, where a clear agreement of the exact and effective dynamics can be observed. In the lower panel, we represent the time-evolution of the coherences, which are damped due to the photon scattering.
More interesting to our purposes are the results displayed in Fig. 17(b), where we have kept the same parameters as above, but set ω L = 10 −3 ω 0 . This guarantees the absence of two-photon Raman transitions. By switching on the additional laser beams, such that Ω 2↓ /ε r = −Ω 1↑ /ε r = 0.05, the non-vanishing differential Rabi frequency in Eq. (A13) should lead to the σ z dipole force (A12). Accordingly, we expect to find damped Rabi oscillations in the coherences for |ψ 0 = (| ↑ + | ↓ )/ √ 2, whereas the populations should only show an exponential damping. This agrees with the results displayed in Fig. 17(b), which also support the accuracy of the effective master equation in this regime.
To recover the action of the lasers on the vibrational degrees of freedom, one should substitute Ω l,s → Ω l,s e ik l ·r in all the expressions above. Besides, in the regime of interest (Ω eff /π)t 0 6 (Ω eff /π)t 0 6 0 −0.9 P ↓ (t) (Ω eff /π)t 0 6 (Ω eff /π)t 0 6 2∆ Ω 1,↓ Ω * 2,↓ − Ω 1,↑ Ω * 2,↑ . Let us note that the running wave of strengthΩ L that only couples to the vibrational excitations, together with the ac-Stark shifts contained in (A7), must be compensated experimentally. Finally, the expression for the effective jump operators is the following Since we are considering that |Ω l,s | ∆, it becomes clear that the effective scattering rates scale as Γ eff = Γ(|Ω l,s |/∆) 2 . the ion trap, or to the additional lasers used for cooling. In particular, we focus on the transverse phonons since they are responsible for the spin-spin interaction. One term in the effect of micromotion can be understood from the interplay between the excess micromotion and those terms in the vibrational Hamiltonian (3) that do not conserve the number of vibrational excitations. By working in the local vibrational basis used in Sec. III, these terms amount to where the couplingsṼ yy i j (t) = V yy i j (t)/(e 2 /l 3 z ) correspond to those of Eq. (4) after taking into account the micromotion In order to neglect ∆H(t), which is responsible of the r.f. heating, we expand Eq. (B5) to leading order in q α 1, and find that all the relevant terms can be neglected under a RWA if the following condition is fulfilled Since the ladder compounds present κ y 1, and Ω rf ω z , the validity of the RWA is easily fulfilled. Thus the leading r.f. heating mechanism would be due to non linearities. This mechanism was analysed numerically in [101,102], where scaling laws for the r.f. heating rates were predicted. In these studies, the heating rates originate from non linearities and thus scale with the initial temperature of the crystal. For the temperatures that apply for our scheme, the r.f. heating would be much smaller than the anomalous heating and thus could be neglected in our analysis.
A different possibility is that of laser heating. Note that the QS requires laser cooling of the transverse phonon modes, although not necessarily to the ground-state. As emphasized in [72,106], depending on the ratio of the r.f. frequency to the decay rate of the cooling transition, either a broadening of the transition or the appearance of multiple micromotion sidebands may occur, which can lead to undesired heating even when the laser frequency is tuned below the atomic resonance. To overcome this effect, one must carefully tune the laser frequency according to the regimes described in [72].
Micromotion contribution to the spin-dependent forces: An important question to address is whether the micromotion modifies the spin-dependent dipole force, as derived in Appendix A. In the interaction picture, the Hamiltonian describing the laser coupling (A2) for the whole ion crystal becomes 2 Ω l,s |r i s i |e ik l ·r i e iξ li cos(Ω rf t) e iδ l,s t + H.c., (B8) after including the micromotion. Here, we have introduced ξ li = q x k l ·r 0 i /2, which represents the ratio of the radial excess micromotion (B3) to the wavelength of the laser radiation.
To carry on with the analysis, we need to specify a particular laser-beam arrangement (see the inset of Fig. 2(b)). We parametrize the laser wavevectors as follows k l = 2π λ sp (cos α l e x + sin α l e y ), where α l determines their angle with respect to the x-axis, and λ sp is the wavelength of the n 2 S 1/2 -n 2 P 3/2 transition (see Fig. 3). Let us recall that in order to control the anisotropy of the spin interactions (19),J eff i j ∝ cos(φ i j ), the corresponding angles must span the range φ i j ∈ [0, 2π]. This implies that the laser-beam arrangement must fulfill α 1 = −α 2 − ∆α, where ∆α |α l |. One particular choice is the following where d is the inter-leg distance (see Fig. 4). With this choice, we find ξ 2i = 0, and ξ 1i = ± π 4 q x , both fulfilling |ξ li | 1. This property will allow us to truncate the following series e iξ li cos(Ω rf t) = ∑ m∈Z i m J m (ξ li )e imΩ rf t , where J m (x) are the Bessel functions of the first kind. By substituting the series (B11) in the laser-ion coupling (B8), we obtain a sum over all possible micromotion sidebands 2 Ω l,s J m (ξ li )|r i s i |e ik l ·r i e i(δ l,s −mΩ rf )t + H.c..
(B12) According to Eq. (A15), we can set ∆ ≈ δ l,s , such that for a sufficiently large detuning Ω rf /2π ≈ 0.1 GHz ∆/2π ≈ 10 GHz, we find that the micromotion sidebands only introduce a resonance for m ≈ 10. However, the contribution of such terms is negligible since ξ li 1, and J m (ξ li ) ∝ (ξ li ) m . As shown in [103], the correction to the leading term should be taken into account by adjusting the Rabi frequency. Let us also note that the remaining non-resonant sidebands can also be neglected in a RWA, since |Ω l,s | ∆. Hence, we conclude that the validity of the spin-dependent dipole force derived in Appendix A is not compromised by the micromotion.
To quantify the accuracy of this argument, we introduce where σ x (t) = Tr |↓ ↑|e −iω 0 t ρ(t) + c.c. represents the coherences. In this expression, the effective evolution under the dipole force σ x (t) eff , as represented in Fig. 17(b), is compared to the time evolution including the micromotion sidebands σ x (t) mic . Accordingly, ε m sets an upper bound to the error of the dipole force (10). In Fig. 18, we represent this error bound as a function of the relative micromotion amplitude ξ 1i . We use the same parameters as in Appendix A, and set the r.f. frequency to Ω rf /ε r = 5 · 10 −3 , which is consistent with the constraint Ω rf ∆. We observe that the micromotion sidebands only contribute with a small error (4-5%) for the regimes of interest ξ 1i = π 4 q x , ξ 2i = 0 (shaded region). Once this has been shown, we must consider the effects of the micromotion on the derivation of the spin-phonon term (12). Following a similar procedure, we obtain the micromotion contributions to the Lamb-Dicke expansion connecting the dipole force (10) to the spin-phonon coupling (12)  (B15) Note that the above expression (B14) includes all the possible resonances between the secular and micromotion sidebands with the laser beatnote ω L = ω 1 − ω 2 . The leading-order term for small Lamb-Dicke parameter η n⊥ 1 occurs for the secular resonance m = 0, and p = 1. This term leads directly to the desired spin-phonon coupling (12) provided that Considering the different orders of magnitude in the problem ω L /2π ≈ ω y /2π ≈10 MHz Ω rf /2π ≈ 0.1 GHz, the leading micromotion resonance would occur for the term Ω L η p n⊥ J m (ξ Li )e i(−pΩ ⊥ n +mΩ rf −ω L )t with m = 1 and p ≈ 9. Since these terms scale as (η n⊥ ) p with η n⊥ 1, they get exponentially suppressed and can be thus neglected. Let us note that the remaining off-resonant sidebands can also be neglected via a RWA in the regime of interest Ω L ω L Ω rf . Therefore, we conclude that the micromotion does not modify the spin-phonon coupling (12).
Micromotion contribution to unwanted transitions: Equation (B8) implicitly assumes that the dynamics due to the Raman-beam configuration can be accurately described by only three levels |↑ i , |↓ i , |r i . While this is always true for Zeeman ions, some special care must be taken for hyperfine ones, where the laser beams may excite some of the remaining states of the ground-state manifold due to the extra resonances introduced by the micromotion. Hence, the laser-ion Hamil-tonian in Eq. (B8) must be supplemented V + ∆V by ∆V (t) = ∑ l,i,a Ω l,a 2 |r i a i |e ik l ·r i e iξ li cos Ω rf t e iδ l,a t + H.c., (B17) which includes all the additional states |a i of the ground-state manifold with energies ε a , and we have introduced the detunings δ l,a = ε r − ε a − ω l . These detunings can be controlled by the Zeeman shifts of the ground-state manifold {|a i } caused by an external magnetic field. Note that the aforementioned unwanted transitions that take the state out of the spin subspace s ∈ {↑, ↓} follow from two-photon processes evolving like Ω 1,a Ω * 2,↑ J m (ξ 1i )J m (ξ 2i ) * e −i(δ a,↑ −ω L −(m−m )Ω rf )t . Therefore, the Zeeman splittings δ a,s = ε a − ε s must be tuned to |Ω 1,a Ω * 2,s | |δ a,s − ω L ± Ω rf |, such that these unwanted transitions become highly offresonant and can be neglected in a RWA. where T is the temperature determining the phonon Gibbs state ρ th = Z −1 e −β ∑ n Ω ⊥ n a † n a n , such that Z = Tr{e −β ∑ n Ω ⊥ n a † n a n } is the partition function, and β = (k B T ) −1 is expressed in terms of the Boltzman constant k B . By considering a separable initial state ρ(0) = |ψ s ψ s | ⊗ ρ th , where |ψ s is a pure spin state, we find the following expression for the relative thermal error This expectation value can be evaluated exactly to yield the following result where the effects of the zero-point motion have been included in σ x i T =0 , and we have introduced the mean phonon numbers for each of the vibrational modesn ⊥ m = a † m a m . It is interesting to note that, as argued for the so-called quantum phase gates [51], the error can be minimized by considering evolution times that are multiples of the detuning of the closest vibrational mode δ ⊥ m * , namely t f = 2πn/δ ⊥ m * , where n ∈ Z. In order to check the validity of our derivation, we have confronted the prediction (C6) to the exact time evolution of the spin-phonon model in Eq. (C1). We use the parameters of Sec. III, namely, the trap frequency ω y /2π = 20 MHz, and the laser parameters η y = 0.1, ω L = 1.1ω y , Ω L = 0.15|δ y |/η y , and e x · k L = 0. In Fig. 19, we represent the exact results for the thermal error ε ex T (dotted lines) obtained by the numerical integration of the Liouville equation dρ/dt = −i[H p +H d , ρ] up to t f = π/8J eff for the Hamiltonian (C1). Note that we truncate the Hilbert space of each vibrational mode ton t = 6 (yellow circles),n t = 7 (green squares), andn t = 8 (blue diamonds). These results are compared to the analytical estimate ε an T (red solid line) in Eq. (C6), showing a remarkable agreement for a sufficiently large truncation of the phonon Hilbert space. We now derive a useful expression for the scaling of the error in terms of experimental parameters, in the limit of |F in | δ ⊥ n . We perform a Taylor expansion of Eq. (C6) for κ y 1, and take into account the orthonormality properties of the normal-mode displacements M ⊥ in . For low-and high-temperatures, we haven ⊥ n ≈ β Ω ⊥ n andn ⊥ n ≈ (β Ω ⊥ n ) −1 respectively, where β = 1/k B T . In both regimes, we find the following upper bound for the error where we have introducedn y = (β ω y ), which corresponds to the mean phonon number for the center-of-mass mode. Finally, in the low-temperature regime, one finds ε T ≤ ε th = 4|Ω L | 2 η 2 y δ 2 yn y (C8) By controlling the experimental parameters, this error term should be minimized for the QS. Note that this linear scaling of the relative error with the mean number of phonons coincides with the results shown in Fig. 19.
Another important question to address is the error caused by heating mechanisms. The heating of a particular vibrational mode may be induced by a variety of factors, such as the combination of stray electric fields and fluctuating trap parameters, elastic collisions with a background gas, fluctuating patch fields in the trap electrodes, or non-linear static electric fields (see [6] for details). Therefore, we shall not focus on a particular microscopic model, but use instead a phenomenological master equation where the coherent dynamics is given by Eq. (C1), and the heating dissipator corresponds to D h (ρ) = ∑ n Γ h a † n ρa n − 1 2 a n a † n ρ − 1 2 ρa n a † n .