Dual trapped-ion quantum simulators: an alternative route towards exotic quantum magnets

We present a route towards the quantum simulation of exotic quantum magnetism in ion traps by exploiting dual relations between different spin models. Our strategy allows one to start from Hamiltonians that can be realized with current technology, while properties of an exotic dual model are inferred from measurements of non-local, string-order-like, operators. The latter can be achieved from fluorescence, or from certain spectroscopic measurements, both of which can be combined with finite-size scaling by controlling the number of ions in the dual quantum simulator. We apply this concept to propose quantum simulators of frustrated quantum magnets, and Ising models with multi-spin interactions. We test the validity of the idea by showing numerically that the predictions of an ideal dual quantum simulator are not qualitatively modified by relevant perturbations that occur naturally in the trapped-ion scenario.


I. INTRODUCTION
Nowadays, the importance of the Ising model in the theory of statistical mechanics can be hardly exaggerated [1]. Nonetheless, early results proved the absence of spontaneous magnetization in the one-dimensional (1D) model and pointed to the lack of phase transitions also in higher dimensions [2], thus stimulating the appearance of other models to explain ferromagnetism [3]. The arguments against the existence of phase transitions in higher dimensions were proven wrong by * Electronic address: bermudez.carballo@gmail.com subsequent efforts [1], such as the exact solution of the Ising model on a 2D square lattice [4], and eventually turned the Ising model into a paradigm of critical phenomena.
Moreover, the relevance of the Ising model goes far beyond the theory of thermal phase transitions. For instance, the partition function of the square-lattice Ising model, and thus its critical phenomena, can be related to the ground state of the 1D Ising model subjected to a transverse magnetic field that plays the role of the temperature [5]. This leads to a new kind of critical effects, i.e. quantum phase transitions [6], where the ground state of a system changes abruptly as a microscopic parameter controlling the quantum fluctuations is varied. Thanks to its exact solution [7], this quantum Ising model also plays a pivotal role in the theory of quantum criticality [6] and, despite its apparent over-simplified appearance, has turned out to be a faithful representation of certain magnetic materials [8].
The Ising model also occupies a privileged position in the theory of emergence in many-body systems [9]. In certain geometries, such as the triangular lattice [10], antiferromagnetic Ising interactions cannot be simultaneously minimized, an effect known as frustration. In the 3D pyrochlore lattice, this magnetic frustration resembles the ordering of protons in water ice [11], and inhibits the development of magnetic order even at very low temperatures. This leads to exotic emergent effects [13] that can be measured in these so-called spin-ice materials [12], constituting a hallmark in the physics of spin liquids. Moreover, if quantum fluctuations are introduced, the wisdom is that even more exotic properties survive at zero temperatures, leading to the elusive quantum spin liquids [14].
From this perspective, it is thus important to understand the interplay of frustration and quantum fluctuations in the Ising model. The standard approach to introduce quantum fluctuations on the vastly degenerate Ising ground state considers additional Heisenberg exchange couplings of tunable strength, which flip pairs of spins and lead to the so-called XXZ model [15]. However, this is not the unique possibility. For instance, one may introduce quantum fluctuations by a transverse field leading to frustrated quantum Ising models, which can also host quantum spin liquids in certain lattices [16]. As a starting point towards this goal, envisaged arXiv:1511.04323v1 [cond-mat.quant-gas] 13  The combination of phonon-mediated spin-spin interactions in linear ion chains with the ability to measure highly non-local operators would allow to exploit a quantum duality transformation to explore exotic Ising models with tunable frustration or multi-spin interactions that can be defined in two-leg ladders with triangular motifs.
in this paper, it is useful to identify setups where the physics of simpler frustrated quantum Ising models can be explored experimentally.
In addition to frustration, interesting many-body effects may also arise as one departs from the standard scenario of pair-wise interactions, and considers multi-particle couplings. For instance, the partition function of a spin-ice Ising model that combines 2-and 4-body interactions in a square lattice [17] can be mapped onto the so-called eight-vertex model [18], which is a paradigm for critical phenomena in ice-type models. In contrast to the standard Ising model with pairwise interactions [4], this Ising model with multi-spin interactions presents some rather exotic properties, e.g. critical exponents vary continuously with the strength of the 4-body interactions, thus connecting different universality classes. In addition to this example, we note that even more exotic phenomena occur when (i) quantum fluctuations or/and (ii) higher dimensions are also considered: (i) the combination of 4-body Ising interactions with a transverse field in the square lattice can lead to phases with topological order [19]. (ii) Ising models on the cubic lattice with 4-body interactions possess a local gauge symmetry that forbids the use of the standard local order parameters [20]. By introducing quantum fluctuations through a transverse field, such models yield a (3+1)dimensional Ising Gauge theory [21,22] through the Euclidean generalization of the quantum-classical mapping [5].
From this perspective, it is important to understand the interplay of multi-spin interactions with quantum fluctuations, and its role in the emergence of exotic quantum magnets. A useful tool to reach this goal, again envisaged in this paper, would be an experimental setup where simple quantum Ising models with multi-spin interactions can be implemented.
In this manuscript, we argue that current experiments of trapped-ion crystals can address the before-mentioned goals: We explicitly demonstrate that relevant instances of frustrated quantum Ising models, as well as exotic quantum Ising magnets with multi-spin interactions, can be implemented in such systems. For our proposal, we exploit a theoretical concept known as quantum dualities [21], which provides a mapping between different spin models. Broadly speaking, a dual map replaces the spin operators at the sites (vertices) of a lattice by analog operators acting on the bonds (links) between the vertices. Accordingly, duality can be used as a tool to switch between local and non-local models, which, as shown in this paper, facilitates the implementation of some Hamiltonians of interest in a trapped-ion setup (see Fig. 1).
The paper is organized as follows. In Sec. II, we describe two possible strategies to simulate the required spin Hamiltonian in linear chains of trapped ions. In Sec. III, we review the notion of quantum duality, highlighting its potential for quantum simulators such as trapped-ion experiments, which allow for highly non-local measurements. This quantum duality is applied to propose quantum simulators of exotic quantum Ising magnets that incorporate the effects of frustration in Sec. IV, or multi-spin interactions in Sec. V. Finally, we present our conclusions and outlook in Sec. VI.

II. TRAPPED-ION LONGITUDINAL XY MODEL
The ever-improving experimental control of certain quantum systems, especially in the realm of atomic, molecular, and optical physics, has reached such a level that it is now possible to devise systems that will behave according to specific quantum many-body models [23]. This synthetic quantum matter, usually known as a quantum simulator, provides an alternative route to address longstanding problems in condensed-matter physics, such as the interplay of frustration and quantum fluctuations mentioned above. Among the different experimental platforms available nowadays, small crystals of atomic ions confined by electromagnetic traps [24] have already proven to be quite flexible quantum simulators [25]. Let us briefly review the possibilities of trapped-ion quantum simulators in the context of frustrated and multi-spin quantum Ising models.
Following the seminal proposal to use the phonons of the ion crystal as mediators of spin-spin interactions [26] (for earlier related proposals see [27]), small-scale 1D quantum Ising magnets have been experimentally realized [28]. The most direct approaches to include frustration would require engineering 2D ion lattices with triangular motifs, either by considering Penning traps [26,29] and micro-fabricated trap arrays in a surface electrode [26,30], or by coping with the excess micro-motion associated with 2D lattices in common rf-traps [31]. A different approach, which has proven very fruitful [32], retains the 1D geometry of the crystals in linear Paul traps while exploiting the long-range character of the Ising interactions to induce frustration via next-to-nearestneighbor couplings. In this work, we shall introduce an alternative method that also permits retaining the 1D lattices, but offers more flexibility by allowing full control of the range of frustration in the simulated spin models.
Regarding multi-spin interactions, there have been propos-als [33] that generalize the above schemes [26] by considering non-linear spin-phonon couplings that can induce 2-and 3body Ising interactions at the expense of getting weaker spinspin couplings with respect to the pair-wise scenario. Hence, synthesizing other multi-spin interactions, such as the mentioned 4-body Ising couplings, will result in even weaker coupling strengths that question the possibility of an experimental observation. A different and very fruitful approach is to build stroboscopic quantum simulators by concatenating wellcontrolled unitaries [34]. For the particular unitaries available in trapped-ion setups, multi-spin interactions can be obtained by introducing simple sequences that involve auxiliary spins [36], as demonstrated in the experiments [36,37]. In this work, we propose an alternative method to implement 4body interactions in an analog fashion, thus not limited by the accumulated error of the concatenated unitaries, but without the above-mentioned limitation on the coupling strengths. The central idea of this work is to exploit the concept of quantum duality [21] to reach the desired quantum simulator of exotic Ising models with frustration, or multi-spin interactions, starting from a different Hamiltonian that is easier to realize in a certain setup. In particular, we will focus on the following physical Hamiltonian which shall be referred to as the longitudinal XY model ( XY), and is the one to be implemented in the ion-trap setup. Let us remark on the two unique features that make this synthetic magnet very different from real magnetic materials: (i) the exact number of spins can be experimentally controlled, (ii) the magnetization, or any other correlation function, can be measured locally. These two features, in combination with the quantum duality, will be crucial for the quantum simulation scheme hereby proposed. Let us now describe the required ingredients to realize the above Hamiltonian (1), and postpone the duality connection to the desired Ising models to the following sections. We start by describing the standard approach in Sec. II A, and then introduce an alternative that might be advantageous in Sec. II B.

A. Spin model from a pair of spin-dependent forces
We consider a collection of laser-cooled atomic ions confined by a linear Paul trap, and forming a 1D Coulomb crystal with equilibrium positions r 0 i = z 0 i e z for i ∈ {1, · · · , N} [38] (see Fig. 1). The small vibrations around these equilibrium positions can be described in terms of three phonon branches, two radial α ∈ {x, y} and one longitudinal α = z, namely H p = ∑ n,α ω n,α a † n,α a n,α , where we have introduced the bosonic operators a † n,α , a n,α that create-annihilate a phonon associated to the normal-mode frequency ω n,α , where n ∈ {1, · · · , N}.
The spin chain is formed by a pair of long-lived electronic levels from the energy-level structure of each particular ion, denoted as |↑ i , |↓ i . This leads to a trivial spin Hamiltonian where we have introduced the transition frequency ω 0 , and the Pauli spin operator σ z i = |↑ i ↑ i | − |↓ i ↓ i |. These spins can be coupled by phonon-mediated interactions, which require a particular form of spin-phonon coupling. One typically considers the so-called spin-dependent dipole forces, which follow from a moving optical lattice, and read as follows where we work in the interaction picture with respect to Eqs.
(2)-(3), and have neglected rapidly oscillating terms within a rotating wave approximation. Here, we have introduced the two remaining Pauli spin operators σ , the zero-point motion of the different vibrational modes ∆r n,α , the dipole light forces F α,β i,n , and their corresponding detunings ∆ n,α,β . We shall assume that a particular light force will only couple to a single vibrational axis, namely F α,β i,n = F α i,n δ α,β . If the trap frequencies are different ω x = ω y = ω z , such that their differences are much bigger than the far-detuned dipole forces |F α,β i,n ∆r n,α | |∆ n,α,β | |ω α − ω β | [26], the phonons act as mere carriers of spin-spin interactions where the phonon-mediated spin-spin couplings have the standard second-order scaling O(F 2 /∆), and thus correspond to virtual processes where a phonon is created and subsequently absorbed by a pair of distant spins in the chain. We are interested in anisotropic XY models (1), and thus shall only consider dipole forces that couple to the two radial branches F x i,n , and F y i,n , such that the effective Hamiltonian (5) only contains XX and YY couplings. In addition, we shall consider a so-called longitudinal field h, which can be implemented by considering the standard carrier transitions [24] with the correct phase. The combination of these ingredients yields the desired XY model of Eq. (1) with the corresponding spin-spin couplings in Eq. (6). As announced above, we shall not exploit the long-range nature of such couplings to engineer frustration [32]. Instead, we will consider the largest possible detunings, which lead to the smallest possible errors in the quantum simulation [26], and also yield a fast-decaying dipolar tail of the anti-ferromagnetic spin-spin interactions We thus obtain a dipolar version of the famous XY model [39], subjected to an additional longitudinal field. B. Spin model from a single spin-dependent force Let us note that, so far, there has not been any experimental realization of a spin model exploiting simultaneously two phonon branches of a trapped-ion crystal. Instead, the dynamics of excitations in the isotropic limit of the XY model has been explored by considering a single Ising model, and thus a single phonon branch to mediate the interactions, subjected to an additional very strong transverse field [40,41]. Unfortunately, this isotropic limit leads to J x i j = J y i j in the above Hamiltonian (1). As will become apparent in the following sections, this would limit its interest as a simulator of exotic quantum magnets under the duality transformation. We now introduce a scheme to modify the implementation of Refs. [40,41] and obtain a fully-tunable XY model in a longitudinal field by exploiting a single phonon branch.
Let us consider the Ising model in a longitudinal field, which arises from considering a spin-dependent dipole force coupled to a single radial branch and a carrier transition Instead of applying a strong transverse field to obtain the isotropic XY model [40,41], let us consider a periodicallymodulated transverse field, which may arise from the crossbeam ac-Stark shift from a pair of co-propagating lasers with different frequencies. This is described by where Ω d is the strength of the transverse field, and ω d its modulation frequency, although we note that the following argument could as well be applied for a transverse detuned carrier by simply exchanging σ z i → σ y i . We start by moving to an interaction picture with respect to the driving, to obtain the following dressed Hamiltonian . By using the Jacobi-Anger expansion of the exponential in terms of Bessel functions, this modulation function can be expressed as We can now simplify the Hamiltonian considerably if we consider that the modulation frequency is much larger than the couplings of the Ising model (8), namely |h|, |J x i j | ω d . By applying a rotatingwave approximation, we find a time-independent Hamiltonian where one observes that some of the coupling strengths are dressed by the zero Bessel function. By moving back to the standard Pauli spin operators, we obtain the desired XY model (1), H I (t) = H XY , with the following parameters Hence, by controlling two ratiosh/J x 0 , and Ω d /ω d , we can access the full phase diagram of the target XY model (1). We note that this scheme might be considered as a spin version of the so-called coherent destruction of tunnelling for electrons in solids under periodic electric fields [42], and is the simplest possible modulation scheme that can lead to interesting effective Hamiltonians. In the context of trapped ions, other driving terms have been exploited to simulate the effects of synthetic gauge fields in the vibrational or spin sectors [43,44].
Before moving to the next section, where we exploit the tool of quantum dualities for quantum simulations, let us summarize the ingredients that we have introduced so far. (i) By loading the ion crystal in a controlled fashion, we can design the number of spins N in the synthetic quantum magnets as desired. (ii) By combining the resonance fluorescence in a cycling transition [24] with global single-spin rotations in the Bloch sphere (i.e. single-qubit gates), we can measure the observables , and so on. (iii) By using carrier transitions and far-detuned spin-dependent dipole forces, either along two different vibrational axes with different trap frequencies, or combining a single force with a periodic modulation, we can obtain a synthetic longitudinal XY model with dipolar couplings (1). These three ingredients will be fundamental to propose a dual quantum simulator in the following section. We shall argue that, even if the dipolar XY model in a longitudinal field does not look more exotic than the standard quantum Ising model, it can lead to very interesting phenomena such as quantum frustration when complemented by quantum dualities and certain local measurements which, although complicated for other setups, are accessible in trapped-ion experiments.

III. DUALITY FOR QUANTUM SIMULATIONS
The notion of duality in physics, which first appeared in connection to the equations of electromagnetism in the absence of sources, has become a far reaching concept in different disciplines. In this work, we shall focus on the use of dualities to understand the properties of many-body lattice models, which started with a seminal work on the self-duality of the Ising model on the square lattice [45]. This property permitted locating the critical temperature of the model prior to its exact solution [4], and turned out to be a powerful concept that can be generalized to other situations [46], including also the quantum-mechanical regime [21,22].

A. Quantum duality transformation
In essence, a duality transformation relates the properties of a particular model in the original lattice Λ with those of a transformed model in the dual lattice Λ * . In 1D, which is the regime of interest to our work, the dual lattice is obtained by placing vertices at the links of the original lattice, and setting bonds between those vertices separated by a site of the original lattice. Therefore, the dual lattice of a 1D chain is simply another chain that has been translated by half a lattice constant. Typically, this transformation is defined for infinite Bravais lattices, and aims at giving a different perspective of the model under study in the so-called thermodynamic limit. However, since we are interested in small-scale synthetic magnets, one needs to consider the effects of boundary conditions, which map an original N-site chain Λ N onto an Once these lattices have been defined, let us introduce the quantum duality on the spin operators [21,22]. The su(2) algebra of the spin operators is preserved by defining the dual spin operators as follows Then, the dual Hamiltonian to the XY model (1) becomes where the sums must be extended to the dual lattice Λ * N+1 , and thus to N + 1 sites. At first sight, this dual mapping does not seem to produce any interesting result, since the dual Hamiltonian seems to be more convoluted than the original one. However, exploiting the fast dipolar decay and some additional properties of the spin-spin couplings (7), this dual Hamiltonian can host a variety of exotic effects. Let us postpone this discussion for the following sections, and address first a question of relevance for the dual quantum simulator.

B. Quantum dual measurements
Usually, quantum dualities are used as theoretical tools that allow to understand certain aspects of a particular lattice model from a different perspective. In this manuscript, we are instead proposing to use them experimentally, which is usually hampered by the non-local character of the dual operators (14). In order to understand the properties of the dual model (15), one typically studies two-point correlators in the dual basis, which become highly-nonlocal correlators in the original lattice, such as Therefore, in order to exploit the quantum duality experimentally and understand properties of the model (15) by performing experiments with the physical Hamiltonian (1), one needs to measure string-order-like parameters, which is usually impossible in most physical setups. However, as argued in the previous sections, these type of observables are exactly the ones that can be obtained via the combination of single-qubit gates and fluorescence in trapped-ion experiments. Apart from correlation functions, also spectral properties give interesting insight into the behavior of a physical system, and important properties such as energy gaps of low-lying excitations are the same for the original model and its dual. Accordingly, spectroscopic measurements in the spirit of the recent experiments [41,47] may provide an alternative way of obtaining relevant information valid for both models.
Before entering into the discussion about the exotic magnetic phases that the Hamiltonian (15) can host, let us also comment on the other highly-unusual property of our synthetic quantum magnets: the possibility of exactly controlling the number of spins N. Critical phenomena and quantum phase transitions take place in the thermodynamic limit (N → ∞). Although this is a mathematical idealization since real systems are always finite, this type of predictions agree considerably well with those based on moderately-large samples, as predicted by finite size scaling (FSS) [48]. Although, strictly speaking, phase transitions can only occur at the thermodynamic limit, FSS can be considered as a window to extract their infinite-size characteristics (e.g. critical points and exponents) by studying the scaling of certain observables with the system size. As occurs with quantum dualities, FSS has been mainly a theoretical tool used in combination with numerical simulations of small systems of increasing size. On the other hand, experimental FSS in real materials would be hampered by difficulties on (i) controlling accurately N, and (ii) distinguishing finite-size effects from other spurious rounding effects (e.g. disorder, impurities, etc.). This situation has recently changed with trapped-ion experiments [28], where the number of Ising spins can be controlled exactly, and where typical disorder and impurity effects are totally absent. Moreover, the possible rounding effects caused by other sources of noise can be experimentally controlled and identified.
In the following two sections, we will combine the duality transformation, FSS, and non-local observables, to show that the mild-looking Hamiltonian (1) can indeed simulate a variety of exotic magnetic phenomena.

IV. DUAL QUANTUM SIMULATORS OF FRUSTRATION
In this section, we want to show that the dual quantum simulator can be used to explore the interplay of frustration and quantum fluctuations. Let us first start with a plausibility argument. One can rearrange the dual Hamiltonian (15) as where is the so-called quantum axial next-nearest-neighbor Ising Hamiltonian (qANNNI), to be defined below. If the XY-model (1) was restricted to nearestneighbour interactions, i. e. if J α i j ∝ δ i, j , the dual mapping between XY and qANNNI model would be perfect, that is, the remainder ∆H = 0 would vanish. However, since our starting point is the XY model with dipolar interactions, ∆H will be non-zero introducing a perturbation to the dual mapping.
The Hamiltonian of the qANNNI model is the one of a transverse Ising model with frustrated nearest-neighbor and next-nearest-neighbor interactions Here, J 1 < 0 is a ferromagnetic nearest-neighbor coupling that corresponds to the longitudinal field of the original Hamiltonian (1). This ferromagnetic interaction competes against an antiferromagnetic J 2 > 0 next-nearest-neighbor coupling, corresponding to the nearest-neighbor XX couplings of the original model, and a transverse magnetic field B, corresponding to the nearest-neighbor YY couplings of the original model. Note that these competing Ising interactions can be described in a two-leg ladder with triangular motifs, which clarifies the origin of the frustration (see Fig. 1). The parameter equivalence under the duality is thus Let us emphasize that the frustration range of the qANNNI model, given by the ratio of next-nearest-to nearest-neighbor couplings J 2 /|J 1 |, corresponds to the ratio of two different couplings J x i,i+1 /h in the original model (1) that arise from totally different sources, i.e. a spin-dependent force and a carrier transition. Accordingly, it is simpler to reach all regimes of physical interest J 2 /|J 1 | ≶ 1/2 by tuning this ratio in the physical model (1), than it would be to reach them by shaping the long-range tail of a bare quantum Ising model [32]. In fact, as J 2 /|J 1 | → 1, it would be no longer justified to neglect the remaining longer-range tail of the bare Ising model, as these additional terms are likely to modify severely the phase diagram. Our method of controlling the range of frustration also seems more straightforward than modifying the crystalline structure in Penning [29] and surface [30] traps, or changing the direction of the lasers/tilting the crystal [31] in common rf-traps.
In addition to the qANNNI Hamiltonian, the long-range nature of the interactions in (1) leads to the following perturbation which contains the dipolar tail of the phonon-mediated interactions starting form the next-nearest-neighbor couplings. Due to the fast dipolar decay (7), the next-nearest-neighbor couplings are almost one order of magnitude smaller than the nearest-neighbor terms in Eq. (18), and the longer-range terms get reduced even further. Hence, one can argue that this perturbation will not modify in any essential manner the phase diagram of the qANNNI model (18). Although the duality between the nearest-neighbor models is already known to bear accurate results, we will show below by numerical calculation that our plausibility argument is correct, and that the additional long-range contributions do not modify the usefulness of the duality mapping for quantum simulation purposes. The model in Eq. (18) is the quantum-mechanical version of the classical ANNNI model [49,50] in a square lattice, as can be proved by applying the quantum-classical mapping [5] in the present case [51]. The classical 2D model, which describes stacked Ising chains with competing ferromagnetic nearest-neighbor and anti-ferromagnetic next-nearestneighbor couplings at finite temperatures, is considered to be a paradigm in the physics of magnetic frustration and commensurate/incommensurate phases (i.e. spatially-modulated arrangements of the magnetic dipoles with a period that is commensurate/incommensurate with the lattice) [52]. In addition to the typical second-order phase transitions describing orderdisorder phenomena in Ising magnets, this model displays different lines separating the commensurate/incommensurate phases, such as Kosterlitz-Thouless [53] and Pokrovsky-Talapov [54] phase transitions, or a disorder line [55] connecting unmodulated/modulated disordered paramagnets. All these thermal phenomena find an analogue in the ground state of the quantum-mechanical model (18), where the transverse field B plays the role of the temperature by introducing quantum fluctuations that drive the different phase transitions.
In the following, we make an exhaustive numerical study to assess the importance of the dipolar terms in the trapped-ion realization, and the possibility of observing such a rich phase diagram with state-of-the-art experimental conditions.

A. Large ground state degeneracy
A hallmark of frustration, sometimes even regarded as its defining property, is the existence of a vastly-degenerate ground state manifold that a arises from the impossibility of minimizing simultaneously the conflicting interactions [13]. In the dual model (18), this frustration arises from the com-petition between ferromagnetic nearest-neighbor and antiferromagnetic next-nearest-neighbor interactions. Triads of sites (i, i + 1, i + 2) define a triangular motif, where the Ising interactions for J 2 /|J 1 | = 1/2 at B = 0 cannot be simultaneously minimized, and thus become frustrated. In fact, the number D qANNNI (N) of degenerate ground states at this frustration point [56] grows exponentially with system size N According to the duality mapping, the nearest-neighbor limit of the XY model in a longitudinal field (1) should show the same exponential growth in the ground state manifold for h = 0, and J x i,i+1 = 2J y i,i+1 . Note, however, that the duality transformation (14) removes the Z 2 -symmetry of the Ising model, and thus halves the number of degenerate states in the original model As discussed above, for open boundary conditions, one also needs to take into account the fact that the duality maps a system of N spins onto a system with N + 1 bond operators. We also note that it is possible to preserve the Z 2 symmetry by modifying the duality transformation such that the transverse field is not applied to all spins [57], but here we stick to the standard duality transformation (14).
We have verified the counting of Eq. (22) for small systems up to N = 14, where the XY interaction is truncated to nearest neighbors. Turning our attention to the dipolar XY model (1), the addition of the long-range tail is especially disturbing around this maximally-frustrated regime, as the vast number of degenerate ground states turns any perturbation, no matter how small, into a highly non-perturbative problem. Nevertheless, the dipolar system still exhibits some traces of the exponentially large degeneracies in the "clean" model. As demonstrated in Fig. 2 by contrasting the 350 lowest eigen-energies of the dipolar model to the energies of the nearest-neighbor model, the dipolar system exhibits a lowenergy manifold separated from higher states by a spectral gap, and the dimension of this manifold is precisely given by Eq. (22).

B. Phase diagram for finite chains
We now explore numerically the phase diagram of the dipolar XY model in a longitudinal field (1), assessing how far it resembles the qANNNI model (18). Let us thus review the features of the qANNNI phase diagram, which are in one-toone correspondence with the 2D classical ANNNI model [52].
Setting B = 0, the qANNNI Hamiltonian becomes classical, and the phase diagram is derived easily: For J 2 /|J 1 | < 0.5, the system is in a ferromagnetic (FM) phase, with a two-fold degenerate ground state having all spins polarized along z. For J 2 /|J 1 | > 0.5, the system is the the anti-phase (AP) with a 4-spin unit cell ↑↑↓↓. For commensurability purposes, one should take a system with a number of spins N divisible by 4. For open boundaries, one obtains a two-fold degenerate ground state (Z 2 spin-inversion symmetry), while in the case of periodic boundaries, translational symmetry yields a fourfold degenerate ground state (Z 2 ×Z 2 symmetry by the Cartesian product of the spin-inversion symmetries at even and odd lattice sites). The critical point between the FM phase and the AP corresponds to the frustration point J 2 /|J 1 | = 0.5 addressed in the previous section, where the model yields an exponentially-large ground state manifold given by (21). If a strong magnetic field B |J 1 |, J 2 is present, the system will exhibit a paramagnetic (PM) phase, with two-point spin correlations that decay exponentially fast with the distance.
Several studies, including exact diagonalization of an effective domain-wall description [58], finite-size scaling [59], or recent matrix-product state (MPS) calculations [60] predict the direct vicinity of the PM and FM phase, with a second-order quantum phase transition along a critical line that bounds the entire FM phase. The FM-PM transition belongs to the same universality class as the standard quantum Ising model [7]. As shown by the density-matrix renormalization group (DMRG) study of Ref. [61], a so-called disorder line, coinciding with the Peschel-Emery line [55], divides the PM phase into two regimes: In direct vicinity to the FM phase, the PM phase is unmodulated, while for larger values of B the exponential decay of correlations is accompanied by a modulation of the correlation function.
For J 2 /|J 1 | > 0.5, a so-called floating phase (FP) separates the AP and the PM phase, as evidenced by the quantum Monte Carlo study in Ref. [62], by the DMRG studies in Ref. [60,63], or by perturbation theory calculations in Ref. [64]. This intermediate phase is characterized by modulated, algebraically decaying spin-spin correlations. The MPS data [60] suggest that the AP-FP commensurate-incommensurate Pokrovsky-Talapov transition is second-order, while the FP-PM transition is believed to be of the Kosterlitz-Thouless type [60,63]. We note that the determination of the phase diagram in this regime is difficult, as important quantities like the energy gap strongly depend on the system size.
Once the properties of the qANNNI model have been discussed, we can start addressing the effects of the long-range tail through the numerical study of the XY model (1).

Parameter region J
In the regime of low frustration, J 2 /|J 1 | < 0.5, the dual qANNNI model behaves similarly to the standard quantum Ising model, with a FM-PM second-order transition. Boundary and/or finite-size effects play a quantitative, but not a qualitative role. The transition between FM-PM has a relatively sharp signature in the spin-spin correlation function at sufficiently large distance. As the PM phase is characterized by exponentially fast correlations, a drop of the correlations marks the boundary. In the original model, a correlation τ z i τ z i+ translates into a non-local -point correlation function σ x i+1 σ x i+2 . . . σ x i+ , which as discussed around Eq. (16) a b c 2 n.n.`XY dip.`XY i Figure (1), and (c) the corresponding correlations τ z i τ z i+ of the associated dual model (18). In all cases, we have applied open boundary conditions, and averaged the spin-spin correlations over all lattice sites. One can see a very nice agreement of these numerical results between the XY model and its qANNNI dual. Note that, in the case of open boundaries, the qANNNI model of N spins corresponds to the dual of a XY model of N − 1 spins, and this is the reason why the spin numbers differ.
As expected, the nearest-neighbor XY model and the qANNNI model behave identically. A comparison of them with the dipolar XY model shows some small differences. For instance, the FM phase is less extended in the presence of dipolar interactions. This is due to the fact that in Eq. (1) the dual FM order is established by the magnetic field term (19), while PM order is due to interactions. Therefore, the dipolar tail, enhancing the interaction, acts in favor of the PM phase and shifts the FM-PM transition towards smaller values of B/|J 1 |. A similar reasoning explains why, along B = 0, the FM phase vanishes for J 2 /|J 1 | < 0.5 in the dipolar model. We also notice that the critical "point" is broadened by dipolar interactions, with a critical region between the FM-AP transition at B = 0.
As a possible approach to determine the FM-PM transition line, we use finite size scaling (FSS) of the mass gap. Therefore, we define the scaled energy gap with z the dynamical critical exponent, in the following taken to be 1. The energies E 1 and E 0 correspond to the first excited state, and to the ground state of the original XY model (1), respectively. At criticality, ∆ N (B) should be independent from the size of the system, and thus all the curves ∆ N (B) for different values N must cross at one common point. This allows us to extract the critical field strength B c for a given value of J 2 . The data obtained in this way, taking into account system sizes up to N = 15, is given by the dots in Figs. 3 and 4. We highlight how the drop of the four-point correlator in Fig. 4 agrees extremely well with the phase boundary obtained via scaling behavior of the mass gap. Let us note that the required energy gap energies can be measured experimentally using the spectroscopic trapped-ion methods put forth in [41,47]. Combined with the unique property of controlling exactly the number of spins, would allow the trapped-ion setups to perform the first FSS in a real experiment.
As an alternative to the energy gap, other observables, such as magnetization and correlations, are equally well suited for performing a FSS, and might be easier to measure in a trapped-ion experiment. Assuming the β = 1/8 exponent of the quantum Ising model, one could scale local magnetization of the ground state as τ z i N 1/8 , which should become independent from N at the phase boundary. In the case of the XY model, we would need to examine the dual correlator to τ z i , that is, ∏ j≤i σ x j . As an example, we have chosen a fixed value J 2 = −0.1J 1 in Fig. 5, and varying the field strength B, we determine its critical value B c by the cut of the curves correspoding to different system sizes. The obtained value, B c = 0.566(10), agrees with the corresponding value obtained via finite-size scaling of the gap, B c = 0.574 (5).
FSS also allows to determine critical exponents by assuming a certain scaling relation [65], namely with ν the critical exponent of the correlation length, B c the critical field strength of the FM-PM transition already calculated, and f (x) a universal scaling function. From Eq. (24), it follows that where ∂ ∆ N ≡ (∂ ∆ N (B)/∂ B)| B=B c . Using Eq. (25), we have determined ν from the data ∆ N with N up to 15, using the a b c n.n.`XY dip.`XY i Figure  B c obtained before. The results are shown in Table I, for different values of J 2 and evaluated using N = 15 and N = 14, with errors assuming an uncertainty of δ B c = 0.005. For all models (qANNNI, nearest-neighbor XY, and dipolar XY), we obtain a critical exponent ν in the range between 0.9 and 1.1, in accordance with the expected Ising universality class ν = 1. Error estimates suggest that the most accurate results are obtained, as could be expected, far off the frustration point J 2 /|J 1 | = 0.5, with δ ν = 0.02. With this, the expected Ising exponent, ν = 1, is located within an interval of up to five times the standard error. We expect that a systematic underestimation of the critical field strength due to the finite size of our samples is responsible for this discrepancy. The best agreement with the Ising value is obtained in the nearestneighbor XY, with values spreading between 0.94 and 1.01. Remarkably, the dipolar version of the XY model produces results which are slightly better than those obtained for the original qANNNI model, see Table I.

Parameter region J
For various reasons, the regime J 2 /|J 1 | > 0.5 is more difficult to capture accurately by studying a small-sized system. First of all, one must be very careful with the number of spins studied, and the choice of the boundary. Since the anti-phase has a 4-spin unit cell, one would take system sizes divisible by 4, such that the classical ANNNI model (B = 0) shows the expected FM-AP transition for both open and periodic boundary conditions at J 2 /|J 1 | = 0.5. However, when we turn to the XY model, the situation is different: If we again choose N divisible by 4, only the system with periodic boundary conditions exhibits a direct FM-AP transition. In contrast, the system with open boundaries, which is the experimentallyrelevant case for the small-scale trapped-ion magnets, shows an intermediate phase, extending from 0.5 < J 2 /|J 1 | < 1.
The occurrence of this intermediate phase can be understood from simple energy considerations in the classical system, i.e. in the XY model without YY interactions. The ferromagnetic case is characterized by a single, fully polarized ground state, whereas the antiphase has Néel order. We assume an intermediate state with ferromagnetic order for N f subsequent spins and Néel order for N a = N − N f spins. The energy of this configuration is given by where we have assumed N a < N. For J 2 /|J 1 | < 0.5, the N a term contributes positively to the energy, and thus the ground state configuration must have N a = 0, that is, the system is in the FM phase. For J 2 /|J 1 | > 0.5, N a should attain its maximal value, which in the formula above is N − 2. The energy for N a = N is simply given by E/|J 1 | = −(J 2 /|J 1 |)(N − 1).
Comparison of the two values shows that the antiphase (i.e. N a = N) becomes favorable only for J 2 /|J 1 | > 1, while for 0.5 < J 2 /|J 1 | < 1, the system has two domains, a ferromagnetic one and a domain in the antiphase. Independent from  the size of the system, the ferromagnetic domain for 0.5 < J 2 /|J 1 | < 1 is of size N f = 2, thus for large systems the intermediate phase becomes indistinguishable from the antiphase.
As we announced before, for a system with periodic boundary conditions, the classical energy is always given by for all possible values of N a . Therefore, the transition from a fully ferromagnetic system to a fully anti-phase system occurs without any intermediate phase.
To solve this disagreement between qANNNI model and XY model, we have to consider the fact that the dual operators are bond operators that reside on the links of the lattice (13). Therefore, the open qANNNI chain of N = 4m spins (m ∈ Z + ) would be the dual model of a Hamiltonian with 4m − 1 spin operators. Accordingly, we should find a much better agreement between the diagram for N = 15 original spins, and N = 16 spins in the dual qANNNI model. This expectation is clearly confirmed by our numerical results (see Figs. 3, 4), which show that the odd number XY model does not present an intermediate phase, but simply the expected FM-AP transition from the dual qANNNI model.
In the quantum case, B = 0, the system is expected to undergo two quantum phase transitions: the AP-FP and the FP-PM transitions. According to an MPS study of the qANNNI model [60], the AP-FP transition is of second order, while the FP-PM is of higher order. Determining the phase boundaries in the thermodynamic limit is difficult, as due to the modulated nature of the FP, the minima of the mass gap strongly depend on the system size N. In principle, as was done in Ref. [61] in the framework of a DMRG study for hundreds of spins, it is still possible to perform finite-size scaling, if one interpolates between different minima of the mass gap ∆ N (B).
Without this costly procedure for extracting the thermodynamic limit, the AP-FP transition at a fixed finite system size can easily be detected by looking at the spin-spin correlators: Sufficiently close to the frustration point, a sudden drop of the 4-point correlations clearly marks a phase boundary in the qANNNI model, the nearest-neighbor XY model, and the dipolar XY model [cf. Fig. 4]. To demonstrate that this drop indeed corresponds to the AP-FP transition, we first note that, in the language of the XY model, the antiphase correlations are characterised by a Néel order of the original spins. Accordingly, the AP spin-spin correlation function is given by c AP (d) := σ x 1 σ x 1+d = (−1) d . We thus introduce the overlap of correlations as For both the nearest-neighbor and the dipolar model with N = 15, the overlap is plotted in Fig. 6  While identifying the anti-phase, and determining its boundary, is possible even for small systems due to the very characteristic correlations of the anti-phase, our method does not provide information about the universality class of the transition. As done in Ref. [60,63], this information can be obtained from the scaling of the entanglement entropy in sufficiently large systems. Such procedure also allows to determine the phase diagram beyond the anti-phase, that is, the position of the higher-order phase transition between floating phase and paramagnetic phase.
Although trapped ions yet do not reach the required system size for identifying the critical behavior at the FP-PM transition, let us nevertheless provide some numerical evidence by studying the scaling of the entanglement entropy in the XY model. Therefore, we first define the entanglement entropy through the nth bond of a system with N spins as S(N, n) = −Tr [1,...,n] ρ [1,...,n] ln ρ [1,...,n] , where ρ [1,...,n] ≡ Tr [n+1,...,N] |Ψ Ψ| the density matrix of the left part of the system (containing spins 1 to n). The system is assumed to be in a pure state |Ψ . Quantum criticality is reflected by a scaling of the entanglement entropy as with c being the central charge of the underlying conformal field theory, and the additive constant being non-universal and depending on n. For non-critical systems the correlation length ξ is finite, and finite-size corrections described by a universal function s(N/ξ ) need to be taken into account. Following Ref. [63], we assume a scaling behavior where the parameter α will be taken as 1. We will also fix n at the center of the chain, that is n = (N + 1)/2, as we will focus on odd N. Eq. (31) implies that S(N) saturates for finite ξ at sufficiently large N, while Eq. (30) implies a logarithmic divergence for critical systems.
Using the DMRG code implemented in the iTensor library [66], we have studied the XY model for system sizes up to N = 327, with up to 200 states kept in a DMRG sweep. Instead of studying the full dipolar model, we have cut the long-range tail of the interactions, taking into account only nearest-and next-to-nearest interactions. For a fixed parameter J 2 /|J 1 |, we calculate S(N) for different magnetic field strengths B. From the behavior of S(N), we are able to determine the critical regions: Up to a first critical field strength B c,1 , the entanglement entropy takes a relatively small value and does not scale with N. This behavior indicates a vanishing correlation length, and agrees with the structure of the anti-phase. For slightly larger values of B, the entanglement entropy still remains constant up to some value N, suggesting that sufficiently small systems still exhibit anti-phase behavior, while larger systems already have floating phase FP behavior. For example, for J 2 /|J 1 | = 1, we get B c,1 ≈ 0.25|J 1 |, but anti-phase behavior persists in systems of size N = 55 for field strengths up to B = 0.26|J 1 |, in systems of size N = 31 up to B = 0.27|J 1 |, and in systems as small as N = 11 even up to B = 0.5|J 1 |.
In the floating point regime, the scaling of S(N) is well described by Eq. (31), compatible with the expected central charge c = 1. Fitting our data to Eq. (31) yields a finite correlation length, plotted in Fig. 7. When B approaches a second critical value B c,2 , the correlation length significantly increases, indicating quantum criticality in the vicinity of the FP-PM transition. For even larger values of B > B c,2 , S(N) again becomes constant, though now in the limit of large N. This behavior is expected for a sufficiently large paramagnetic system, characterized by a small but finite correlation length.
Qualitatively, the described behavior is found for different J 2 /|J 1 |. We have focussed on two values J 2 /|J 1 | = 0.8 and J 2 /|J 1 | = 1. In the former case, we obtain B c,1 = 0.10|J 1 | and B c,2 = 0.45|J 1 |. In the latter case, we obtain B c,1 = 0.25|J 1 | and B c,2 = 0.60|J 1 |. Remarkably, the extent of the floating phase along B, that is B c,2 − B c,1 , is the same for both parameters J 2 , in agreement with earlier findings for the qANNNI model [63].
Finally let us turn our attention to the nearest-neighbor XY model, which we have studied for J 2 /|J 1 | = 1. We find the same behavior as before for the system with nearest and next-to-nearest interactions, but the critical field strengths B c,1 and B c,2 are shifted towards larger values (0.52 and 0.81 for J 2 /|J 1 | = 1). The extent of the floating phase is slightly smaller than in the system with nearest-neighbor and nextnearest-neighbor interactions.
As a summary, our results suggest that, sufficiently far from the frustration point, the dipolar XY model has no qualitative changes with respect to the qANNNI phase diagram. We find clear evidence for the paramagnetic to ferromagnetic phase transition, as well as for the transition from the anti-phase to the floating phase, even in small-sized systems. DMRG calculations for larger systems demonstrate that also the phase transition from the floating phase to the paramagnetic phase persists in the presence of interactions beyond nearest neighbors.

V. DUAL QUANTUM SIMULATORS OF MULTI-SPIN INTERACTIONS
In this section, we want to show that the dual quantum simulator (15) can be used to explore the phase diagram of quantum Ising models with multi-spin interactions.
In order to achieve the desired Hamiltonian, we shall focus on a particular trapped-ion realization of the above scheme that uses two different hyperfine levels within the atomic ground state manifold to encode the pseudo-spin. In this case, the spin-dependent dipole forces (4) that couple to both radial phonon branches, F α i,n with α ∈ {x, y}, arise from a combination of two-photon Raman transitions, each of which requires a pair of laser beams tuned far from the resonance of a dipole- allowed transition to an excited state [24]. If the effective wave vector of the interfering laser beams, ∆k k k α = k k k α L,1 − k k k α L,1 , has a component along the ion-trap axis, e z · ∆k k k α = ∆k α z = 0, and the two-photon Raman transitions lie far away from the axial sidebands, one recovers the desired XY model (1), but with modified spin-spin interactions (7), namely We thus obtain periodically-modulated spin-spin interactions that still decay with a dipolar power law. By controlling the ion lattice spacing, or the direction of the laser beams, one can set ∆k x z z 0 i − z 0 i+1 = π(2m+1) 2 for some integer m ∈ Z + , such that the nearest-neighbor spin-spin couplings are inhibited J x i,i+1 ≈ 0. The fulfilment of this condition requires an homogeneous lattice spacing in the ion crystal, which can be achieved by using micro-fabricated surface traps [30] or by segmenting the axial electrodes in linear Paul traps [67].
Under these conditions, we can rearrange the dual Hamiltonian (15) Here, the Ising interactionJ 2 < 0 is a ferromagnetic nearestneighbor coupling that corresponds to the longitudinal field of the original Hamiltonian (1),J 4 > 0 is an anti-ferromagnetic 4-spin Ising coupling that corresponds to the next-nearestneighbor XX couplings of the original model, andB is a transverse magnetic field that corresponds to the nearest-neighbor YY couplings of the original model. This parameter equivalence under the duality transformation is Let us emphasize that the competition of the Ising interactions, characterized byJ 4 /|J 2 |, corresponds to the ratio J x i,i+2 /h in the original model, which can be tuned to any particular value. The same occurs for the ratioB/|J 2 |, which corresponds to J y i,i+1 /h, such that the whole phase diagram of the qIMS can be addressed with the dual quantum simulator.
In addition to the qIMS Hamiltonian, the dipolar tail of the interactions leads to the following perturbatioñ where some of the couplings also vanish due to the periodic modulation (32), namely J x i j ≈ 0 for all j = i + (2m + 1) with m ∈ Z + . Due to the fast dipolar decay (32), one can argue once more that this perturbation will not modify in any essential manner the phase diagram of the qIMS model (33), such that the dual quantum simulator (15) also gives access to the physics of these exotic quantum magnets.
In the absence of 4-body couplingsJ 4 = 0, this model (33) corresponds to the standard quantum Ising model [7], displaying a second-order phase transition atJ 2 =B in the Ising universality class [4]. In the absence of pair-wise couplings J 2 = 0, the model corresponds to the quantum-mechanical version of a classical Ising model in a square lattice [68] describing stacked Ising chains with 4-body couplings at finite temperatures, as can be proved by applying the quantum-classical mapping [5] in the present case [69,70]. The model has a Z 2 × Z 2 × Z 2 symmetry (i.e. cartesian product of spin inversions for the 3 different spin pairs in each 4-site partition of the lattice), which leads to an eight-fold degenerate ground state forB = 0, and a first-order phase transition atJ 4 =B in the universality class of the q = 8 Potts model [71]. As occurred for the qANNNI model (18), one thus expects that the competition of both Ising interactions, and their interplay with the quantum fluctuations brought by the transverse field, must lead to very interesting critical phenomena. In contrast to the frustrated Hamiltonian, this quantum multi-spin Ising model has not been studied in such detail. FSS studies with up to N = 16 spins already point towards very interesting critical phenomena: The nature of the quantum phase transitions as a function ofB changes from first to second order aroundJ 4 /|J 2 | = 1/2 [72]. However, the limited finite sizes did not allow for accurate studies close to the interesting pointJ 4 /|J 2 | = 1/2, which classicallyB = 0 leads to a largely-degenerate ground state that can give rise to a variety of phases upon switching the transverse magnetic field. We believe that a more careful analysis of this region would be very interesting, and we hope that this manuscript will stimulate future work in that direction.
From the experience gained by the exhaustive numerical study of the previous section, we conjecture that the quantum simulator of the quantum Ising model with multi-site interactions (33) will not be compromised by the additional dipolar tail (35), such that the above interesting region can be characterized experimentally using the accessible non-local order parameters (16), or the FSS of other spectroscopic observables, characteristic of current ion-trap experiments.

VI. CONCLUSIONS AND OUTLOOK
In this manuscript, we have presented an alternative route to build quantum simulators of exotic magnetism by exploiting the notion of quantum dualities. In certain situations, such as the ones discussed in this work, such duality transformations allow one to explore interesting models that involve frustration or multi-spin interactions, and their interplay with quantum fluctuations, by focusing on different spin models that are simpler to implement in a particular experimental platform. This dual approach has one important caveat: measurements of highly non-local observables should be feasible in the particular experiment. This makes trapped-ion setups ideally suited for this playground of dual quantum simulation.
We have focused on two particular dual quantum simulators which can be implemented using state-of-the-art technology in linear ion crystals, and which allow one to explore paradigmatic models of frustration: the quantum axial nearest-neighbor Ising model, and the quantum Ising model with competing 2-and 4-body interactions. In the former case, we have made a careful numerical study to prove the validity of our scheme, which takes into account relevant perturbations that occur naturally in the trapped-ion scenario.
Although we have focused on these two particular examples for ion chains, the notion of dual quantum simulators can be certainly applied to other models, such as the quantum spin liquids, topologically-ordered phases, or Ising lattice gauge theories mentioned in the introduction. Provided that quantum magnets are finally synthesized in Penning or surface traps, this quantum-duality approach will also become relevant for other lattice geometries, which may define an alternative route to these exotic quantum many-body phenomena.
Finally, let us mention that in the recent years another platform for long range spin models has been proposed, namely ultracold atoms in nanostructures, or more precisely ultracold atoms trapped in a vicinity of tapered fibers and optical crystals (band gap materials). The experimental progress in coupling of ultracold atomic gases to nanophotonic waveguides is very rapid (for a recent review cf. [73]). The ideas and proposals concerning realization of long range spin models were developed for instance in Refs. [74][75][76], and it would be interesting to explore the possibilities of quantum dualities in this context.