Trapped-ion quantum simulation of tunable-range Heisenberg chains

Quantum-optical techniques allow for generating controllable spin-spin interactions between ions, making trapped ions an ideal quantum simulator of Heisenberg chains. A single parameter, the detuning of the Raman coupling, allows to switch between ferromagnetic and antiferromagnetic chains, and to modify the range of the interactions. On the antiferromagnetic side, the system can be tuned from an extreme long-range limit, in which any pair of ions interacts with almost equal strength, to interactions with a $1/r^3$ decay. By exact diagonalization, we study how a system of up to 20 ions behaves upon tuning the interactions. We find that it undergoes a transition from a dimerized state with extremely short-ranged correlations towards a state with quasi long-range order, that is, algebraically decaying correlations. On the ferromagnetic side of the system, we demonstrate the feasibility of witnessing non-locality of quantum correlations.


I. INTRODUCTION
A paradigm system of quantum mechanics which may exhibit intriguing quantum properties like entanglement and non-locality are two spins. By increasing the number of spins, more complex behavior may emerge. In fact, a large variety of condensed matter phenomena, ranging from metal-insulator transition to superfluidity or superconductivity, are successfully described by mapping the relevant low-energy Hilbert space onto a spin model [1]. Moreover, spin models may describe spin-liquid phases which exhibit topological order. In recent years, technological progress in manipulating atoms on the quantum level has allowed to explicitly engineer spin models [2]. This has opened the opportunity for testing the foundations of quantum mechanics, and simulating complex many-body behavior.
A very promising quantum simulator are trapped ions. They can be prepared in such a way that their dynamics is mainly restricted to some internal states of the ions, while the external motion is cooled down to only a few phonons. The internal states then represent a (pseudo)spin, and phonon-mediated spin-spin interactions can be implemented [3,4]. This has already led to the experimental realization of SU(2) Ising models in one and two spatial dimensions [5][6][7][8], and, very recently, to the experimental study of entanglement dynamics in Ising and XY chains [9,10]. The implementation of more complicated spin models has been suggested, e.g. of Heisenberg and XY models [4,11], or models with higher spin [12]. Furthermore, the tunability of the phonon-mediated interaction allows to study models with long-range interactions. Roughly, interactions with 1/r α decay have been engineered for 0 < ∼ α < ∼ 3 [8]. This flexibility suggests an implementation of tunablerange spin models in trapped ions. While both theoretical and experimental literature so far has focussed on Ising-or XY-type quantum simulations [5][6][7][8][9][10][11]13], here we consider a trapped-ion implementation of the Heisen- berg model. In particular, by assuming an experimental setup as sketched in Fig. 1 (a,b), we study the influence of a single control parameter on the quantum simulation, the detuning of the Raman coupling. As shown in Fig. 1 (c,d), this parameter controls the range of the interactions. Modifying it may bring our quantum simulation close to different antiferromagnetic variants of the Heisenberg model: the Haldane-Shastry model [14][15][16], the Majumdar-Ghosh model [17], or the Lipkin-Meshkov-Glick model [18]. While the Haldane-Shastry model describes antiferromagnetic order with algebraic decay of correlations, the Majumdar-Ghosh model provides a parent Hamiltonian for a fully dimerized ground state, that is a ground state with extremely short-range correlations. In the Lipkin-Meshkov-Glick model, spinspin interactions are independent from distance, and the system is thus described in terms of global spin operators.
In this paper, we first provide an overview of the different models in Section II. We then discuss the ionic setup in Section III. In Section IV, we show that by increasing the range of interactions, a transition from a Haldane-Shastry-like quasi-long-range order to a Majumdar-Ghosh-like dimerized order can be observed in the ionic system. In Section V, we consider the Lipkin-Meshkov-Glick limit on the ferromagnetic side, which has been suggested for witnessing non-locality of quantum correlations [19]. Finally, we provide a summary and outlook in Section VI. Moreover, our paper contains two appendices: In appendix A, we provide details to the ionic setup, and in appendix B, we discuss the relation between the ionic model and a spin model.

II. THE MODELS
The Hamiltonian of a Heisenberg chain of N spin-1/2 particles generally reads where σ (i) α denotes a Pauli matrix for the spin at position i. In this paper we focus on the antiferromagnetic side of this model, that is, the model with interaction strengths J (i,j) > 0.
The functional behavior of J (i,j) as |i−j| may crucially influence the physics of the model. Let us in the following discuss the different cases which are important for our application.
In this case, interactions decay quadratically with the distance d ij between the two spins, J ij . For an analytic description of the model, it is convenient to impose periodic boundary conditions, and consider spins arranged on the unit circle, that is, with positions z k = exp 2πi N k . The ground state of the model can then exactly be obtained from a Gutzwiller ansatz. For a chain with even number of spins, it reads with M = N/2, S + zi = |↑ ↓| zi a raising operator of the spin at position z i , and the coefficients of each Fock state given by the wave function This function has a remarkable similarity to the Laughlin wave function for two-dimensional systems in the fractional quantum Hall regime [20]. The similarities can be extended to the excited states of the model, which are spinons with certain anyonic properties: Obtained as a superposition of spin flips, excitations live in an integer-spin Hilbert space, but, only occuring pairwise, each spinon carries half-integer spin. In that sense, the spinon represents a quasiparticle with a fractional quantum number. Also, Haldane's generalization of the Pauli principle [21] allows to associate fractional quantum-statistical behavior to the spinons by noticing that each spinon pair reduces the number of available single-particle states by 1, in contrast to fermions which would reduce the number of states by 1 per particle, or bosons where the number of available states would not be affected by the presence of particles. As a fingerprint of Haldane-Shastry-like behavior, one can consider the spin correlations. It has been shown analytically that the model supports correlations with a power-law decay [22] With this criterion, the Haldane-Shastry model, despite the long-range interactions, supports the same quantum phase as the Heisenberg chain with nearest-neighbor interactions, also characterized by quasi long-range spin order. In the next subsection, we will introduce a model which is in contrast to this behavior.
the total state is a product over singlet bonds: It is obvious that in such dimer state, any spin is fully correlated to one nearest neighbor, but fully uncorrelated with the other spins. In other words, no long-range spin correlations exist. Note that for a system with open boundary conditions, the spin at position 1 and the spin at position N do not interact, and therefore the ground state is uniquely given by |Ψ − .
By strengthening interactions between distant spins, one approaches the Lipkin-Meshkov-Glick model [18] in which interactions are spatially independent, that is, J (i,j) = J. With this, the Hamiltonian is rewritten . For antiferromagnetic interactions, any singlet state is thus a ground state, leading to a huge degeneracy. The number of singlet states formed by N spin-1/2 particles is given by the Catalan number However, this degeneracy is lifted by any small spatial dependence of the interactions. As we will see below, the ionic systems approaches the Lipkin-Meshkov-Glick limit when the Raman coupling is close to resonance with the center-of-mass phononic mode. By going through the resonance, one is able to swap the sign of the interactions, thus both the antiferromagnetic and the ferromagnetic Lipkin-Meshkov-Glick model can be realized. On the ferromagnetic side, the ground states are the Dicke states, that is, symmetric superpositions of all states with a fixed spin polarization, that is with a fixed number N ↑ (N ↓ ) of ↑-(↓-)spins. The unnormalized Dicke states read Apparently, the ground state degeneracy is (N + 1)fold, and, as each Dicke state has a different spin polarization, it can be lifted by a polarizing field term z . While such term will make the fully polarized state |↓ . . . ↓ the unique ground state, one can also obtain the polarized Dicke state D N/2,N/2 as the unique ground state by reducing the Heisenberg XXX interactions to XX interactions. The Hamiltonian then reads with D N/2,N/2 the unique ground state for J < 0 and h = 0. .

III. THE IONIC SYSTEM
Tunable-range Ising models have already been implemnted in trapped ions [5,6]. In these experiments, twolevels ions are confined to a line by a Paul trap, and their motional state is cooled to only a few phonons. Strong spin-spin interactions are then achieved by applying state-dependent forces on the ions. In Ref. [5], Ising-type interactions of the form σ z are generated, whereas Ref. [6] describes the production of interactions of the form σ y . As reviewed in Ref. [24], both approaches have basically the same footing, and it is possible to combine them. Instead of an Ising coupling, one then obtains, in the first place, an XY Z-model. By making all interactions equal, one gets the Heisenberg model.
An important difference between the ionic system and the ideal models discussed above are the boundary conditions: For the most feasible experimental implementation, they are open, while periodic boundary conditions are convenient for a theoretical description. In general, the effect of boundary conditions is minimized by scaling up the system. For systems of up to N = 20 ions, we will in the following discuss how close the connection to the Haldane-Shastry model and the Majumdar-Ghosh model can be made by tuning the range of the interaction.
To this goal, let us first briefly review how the desired spin-spin interactions can be engineered. For each coupling, i<j σ α , a pair of Raman lases is set up, as depicted in Fig. 1(a,b). Making several assumptions, which are sketched in the appendix A and detailed in Ref. [24], the Hamiltonian for the interaction of the ions with the photons is given by for α = x, y, and Here,hω 0 is the energy difference between the two internal levels,â m is the annihilation operator for phonons denoted by m and with frequency ω m . The Rabi frequencies of the couplings are denoted by Ω α , and ω α are the frequencies of the fields. Furthermore, the strength of each coupling depends on the Lamb-Dicke paremeters η m , which are explicitly defined in the appendix A. For all couplings α, the wave vector difference of the photons, δk α , is assumed to be transverse to the ion chain, so the sum over m reduces to a sum over N transverse modes. It has been shown in Ref. [6], for a system with a single coupling term h α , that the time evolution can be made identical to the one of a spin system with the Hamiltonian where the spin-spin interaction strength is given by Here,ω α ≡ ω α − ω 0 for α = x, y, whereasω α ≡ ω α for α = z. Some details of the derivation of this formula are provided in the appendix B. In particular, we generalize to the case of more than one coupling, and show that the couplings can be chosen such that they do not interfere. The Hamiltonian then is effectively given by with all J (i,j) α given by Eq. (12). These interactions can be tuned by varying the frequency of the Raman laser. Choosing it close to the frequency of the centerof-mass mode, the strength of the induced interactions will barely depend on the particles' position, J (i,j) ≈ J. The sign of J can be made positive (negative) by tuning above (below) the center-of-mass frequency. Close to resonance, the system is similar to the Lipkin-Meshkov-Glick model. Due to the strong nearest-neighbor and next-nearest-neighbor interactions this limit is somewhat similar to the dimerized J 1 − J 2 -model.
By increasing the detuning from the center-of-mass mode, J (i,j) is made spatially dependent, and in the limit of large detuning, a 1/r 3 decay can be achieved. For intermediate values of the detuning, the interactions may approximate a quadratic decay for sufficiently small r, as shown in Fig. 1 (c,d). The most dominant interactions then agree quantitatively well with the interactions of the Haldane-Shastry model.
Realistically, interaction strengths of the order of kHz can be achieved [7] at a sufficiently small amount of errors. This is enough to keep the time scales of the simulation faster than decoherences from heating or imperfections. The scalability of such quantum simulation has been discussed in Refs. [7,25].
Further techniques could be applied in order to achieve a better agreement of the ion setup with a particular model, e.g. by individually addressing of the ions [26]. However, here we restrict ourselves to the simplest implementation which already exhibits rich physics. in our study is the detuning δ ≡ω α − ω COM , where ω COM refers to the frequency of the center-of-mass mode, that is, the transverse mode of largest energy. For convenience, we have chosen 171 Yb + ions, with equilibirum distance of 10µm at a radial trap frequency ω trap = 2π×5MHz. As shown in Fig. 1(c,d), for very small detuning, δ = 1kHz, the system is close to the Lipkin-Meshkov-Glick limit: Any pair of ions interacts with almost the same strength. For large detuning, δ > ∼ 1000kHz, the interactions decay with 1/r 3 . For intermediate values, the interactions between near neighbors becomes similar to the Haldane-Shastry interactions.
In Fig. 2(a), we plot the overlap of the ground state of 16 ions with the ground state manifold of the Majumdar-Ghosh model, and the Haldane-Shastry model as a function of the detuning. In the limit of small detuning, the ground state is almost fully dimerized. Increasing the detuning, the ground state of the Haldane-Shastry model becomes more relevant. One has to note, however, that (for 16 particles) this state is far from being orthogonal to the dimerized subspace. In fact, it has itself an overlap of 0.73 with the dimerized manifold.
A sharp criterion for the transition from a dimerized state to a long-range ordered state can be inferred from Fig. 2(b), where the energy of the lowest singlet and the lowest triplet excitation is plotted. For the J 1 −J 2 model, it is known that the crossing of these energies mark with high precision, even in small systems, the dimerization transition [27]. In our case, the triplet becomes the lowlying excitation for δ > 10kHz. In this context, we also note that, at δ = 1kHz, in total 118 singlet states have lower energy than the triplet state. Still, this is far from C 16 = 1430, the number of singlet states providing the ground state manifold of the antiferromagnetic Lipkin-Meshkov-Glick model. But in contrast, at δ = 5kHz, the lowest triplet state already reaches the 7th position in the energy spectrum.

A. Entanglement entropy.
A dimer is a pair of maximally entangled spins. A measure which localizes the entanglement within a system, and which is thus able to identify dimerization, is the entanglement entropy. This quantity considers bipartitions of the system, and measures the entanglement between the two subsystems. In our case the spins form an onedimensional array, and it is thus natural to consider the bipartitions A = {1, . . . , ℓ} and B = {ℓ + 1, . . . , N }. The entanglement entropy is then a function of ℓ, defined as whereρ(ℓ) is the density matrix of the subsystem A.
The nearest-neighbor dimerized states, |Ψ ± of Eq. (5), are characterized by a strongly alternating behavior of the entanglement entropy, as shown in Fig 3: A cut through every second bond will yield strong entanglement, S(ℓ) = 1, as the spins ℓ and ℓ + 1 are dimerized, whereas on the other bonds no quantum information is shared, S(ℓ) = 0. Such alternations, however, are not present for the Majumdar-Ghosh dimer with periodic boundary conditions [28]. In that case, both dimer states |Ψ + and |Ψ − , will equally contribute to the ground state manifold, and thus completely wash out the pattern. It has been discussed in Ref. [29] for a broader class of spin models that alternating behavior is characteristic for systems with open boundary conditions, though typically with a much smaller amplitude than in the case of the Majumdar-Ghosh chain.
In Fig. 3, we plot the entanglement entropy of different systems with N = 18 spins: the ionic system with detuning δ = 1kHz and δ = 100kHz, an ideal Heisenberg chain with nearest neighbor interactions, and the Majumdar-Ghosh chain. In all cases, we have applied open boundary conditions, and accordingly we find alternating behavior of the entanglement entropy. As expected, these alternations are strongest for the Majumdar-Ghosh chain, and weakest for the nearest-neighbor Heisenberg model. The entanglement entropy of ionic systems with δ = 1kHz behaves very similar to the entanglement entropy of the Majumdar-Ghosh model, proving the dimerized nature of the phase. At δ = 100kHz, the curve alternates less and comes closer to the entanglement entropy of the nearestneighbor Heisenberg model. This shows that the tendency of nearest neighbors to form singlets is still present, but also more remote spins become entangled.

B. Spin-spin correlations.
An experimentally accessible quantity which nicely displays the different entanglement properties are spin correlations. While correlations will, in principle, depend on the position of the spins, we define an average which depends only on the distance d between the spins: To some extent this restores periodic boundary conditions. Furthermore, we perform a finite-size scaling (taking into account all even system sizes from N = 10 to N = 20). The results, for δ = 100kHz, are shown in Fig. 4. The sign of the correlations alternates with odd/even d. The correlations compare well with the correlations of the Haldane-Shastry model, despite the different boundary conditions. The decay is slightly too fast, but through finite-size scaling a slightly better agreement is obtained. One should also note that for large d, comparable to the system size, edge effects are not averaged out by the definition of Eq. (15). The plot in Fig. 4(b) demonstrates that the decay can be modeled by a power-law: C(d) ∝ (−1) d r −α , with α = 1.15 ± 0.05. On the dimerized side, we find correlations as shown in Fig. 5. Since C(d) as defined in Eq. (15) would average out the large dimer correlations with the essentially uncorrelated bonds between two dimers, we have defined where i max = (N − d)/2 for d even, and i max = (N − d + 1)/2 for d odd. With such definition, we find that C odd (d) takes its maximum value −1/4 for d = 1, demonstrating the strong anticorrelations between every second nearestneighbor pair. For larger distances, the value rapidly descreases. The sign alternates for odd/even d. As shown in Fig. 5, the decay takes place exponentially.

V. NONLOCALITY WITNESS
In the previous section, we have discussed one intriguing aspect of quantum mechanics which can be studied in the ionic system: entanglement. A strongly related and particular striking feature of entanglement is the nonlocality of correlations [30]. Nonlocality can be defined as the impossibility of classically simulating the outcome of a local measurement while a second, remote measurement is performed unless some information obtained during this measurement is shared. As first shown by J. S. Bell, nonlocality is witnessed by the violation of some inequalities which local correlations have to fulfill [31]. In general, unfortunately, Bell inequalities depend on various high-order correlation functions, and it is thus difficult to witness nonlocality in a system. In Ref. [19], however, it has been proposed to detect nonlocality by measurement of two-body correlators only. As a specific example, a Bell inequality has been given which is violated by the spin-polarized Dicke state D N/2,N/2 . It reads Here, are the correlations of two measurementsm taken at sites k and l. On each site, one may takem 0 ≡ σ z andm 1 ≡ cos θσ z + sin θσ x . If a θ for which inequality (17) is violated exists, the system must have non-local correlations.
While inequality (17) turns out to be unable to detect nonlocality on the antiferromagnetic side of the ionic setup, a chance for violation exists on its ferromagnetic side. For sufficiently small detuning, the sytem is in the Lipkin-Meshkov-Glick limit of spatially independent interactions, and for J < 0, states with maximal spin are ground state of the Hamiltonian H ≈ J(S 2 − 3 4 N ). The ground state manifold is thus given by all N + 1 Dicke states. To lift the degeneracy, one could apply a magnetic field: A field in z-direction yields either D 0,N or D N,0 as the unique ground state. The state D N/2,N/2 becomes unique ground state if a staggered magnetic field along z-direction is applied. Alternatively, even without applying a symmetry-breaking field, the degeneracy is lifted in favor of D N/2,N/2 by switching off the interactions J (i,j) z → 0. The system is then govered by H LMG given in Eq. (8).
In these cases, non-locality can be proven via inequality (17), as the spin-spin correlation of D N/2,N/2 are given by With this, inequality (17) is violated. To assert the feasibility of such non-locality witness in the ionic system, we have calculated the ground state of a system of 14 ions, with J given by Eq. (12) . Indeed, the system is found in the Dicke state D N/2,N/2 with fidelity > 0.99 for |δ| < 1kHz. By increasing the absolute value of the negative detuning, we find that, while S zz and S zx remain constant, S xx decreases. Thus, for S xx below a critical value S crit xx = N (N − 2)/2, inequality (17) will not be violated anymore. For N = 14, we numerically found S xx = 89.6 > S crit xx at δ = −2.8kHz, while S xx = 73.3 < S crit xx at δ = −3.0kHz. In the same parameter range, the overlap of the ground state with the Dicke state drops suddenly: While it remains above 0.95 up to |δ| = 2.5kHz, it reaches zero for |δ| = 3.2kHz.
In this regime of relatively large detuning, other phonons than the center-of-mass mode strongly affect the spin-spin interactions. This gives rise t very peculiar in-teraction patterns consisting of attractively and repulsively interacting pairs. It would, however, be striking if non-locality of correlations could also be witnessed in a system where interactions are short-ranged. Interestingly, we find that the short-range ferromagnetic Heisenberg chain behaves exactly the same way as the longrange chain does: N + 1 Dicke states form a ground state manifold in which degeneracies can be lifted by magnetic fields, in particular, in favor of D N/2,N/2 by a staggered magnetic field. Certainly, the ionic setup discussed here does not immediately allow for implementation of the ferromagnetic short-ranged Heisenberg model due to the mentioned mixing of different phonon modes, but it could be achieved by switching from transverse to longitudinal phonons as transmitter of interactions. In that case, the center-of-mass mode is lowest in energy, and the negative detuning will not interfere with other modes. Moreover, additional control could be implemented through additional Raman couplings [26]. Alternatively, also atoms in optical waveguides or phononic crystals allow for implementing spin models with controllable interactions, and could be tuned into a ferromagnetic short-range regime [32].

VI. SUMMARY AND OUTLOOK
We have theoretically studied a quantum simulation of the Heisenberg model which is feasible with trapped ions. In particular, we have investigated the influence of a single control parameter, the detuning, on the simulation. This parameter allows to tune the range of the interactions, and we have shown that this can trigger a dimerization transition. To demonstrate dimerization, we have calculated the entanglement entropy and two-spin correlation functions in systems of up to N = 20 ions. The latter can readily be measured with fluorence measurement techniques. On the ferromagnetic side, we find a parameter window in which measurement of two-spin correlation functions is able to witness non-locality. In summary, the trapped-ion quanatum simulation is able to test basic foundations of quantum mechanics, and to study complicated, long-ranged models.
In this context, it could be particularly interesting to measure also the dynamical structure factor in scattering experiments. If the number of particles is sufficiently large, this should allow for identifying the spinons in the excited states. With this, one could demonstrate that a single spin flip consists of two spinons, and thus, that the elementary excitation of the system, a single spinon, carries spin-1/2, a fractional quantum number. In that sense, the ion chain could be used to prove the existence of anyonic quasiparticle in one spatial dimension. and h α describes the interactions with the photons, Here,hω 0 is the energy difference between the two internal levels, σ z is a Pauli matrix. The phonons in the system are counted byn m =â † mâ m , for each mode denoted by m. The different Raman couplings, denoted by the index α, give rise to Rabi frequencies Ω α . The fields have wave vector k α , frequency ω α , and a phase ϕ α . The position of the ions is denoted by r (i) , and the action of the field on the internal state of each ion is expressed by τ (i) α . We choose the polarization of the couplings such that τ (i) Several assumptions on the Hamiltonian parameters can be made to simplify the expression. First of all, in the Lamb-Dick regime we have k α · r (i) ≪ 1, and we therefore approximate e i(kα·r (i) ) ≈ 1 + ik α · r (i) = 1 + m η (i) mα (â m +â m ). In the latter step, we have replaced the ions' coordinates r (i) by the normal modes around their equilibrium position. The scalar product k α · r (i) is then rewritten in terms of Lamb-Dicke parameters η mα . These parameters measure the strength of the couplings between the photons with wave vector k α and the phononic mode m.
To obtain the values of the Lamb-Dicke parameters, we have to calculate the normal modes in direction of k α , which may be either transverse to the ion chain in x-or y-direction, or longitudinal along the z-axis. We obtain the normal modes, denoted by the N × N matrix M α m,i , by diagonalizing the vibrational Hamiltonian K α , Here, we label the phononic modes by m ∈ 1, . . . , N and an additional label α ∈ {x, y, z} specifying the direction where c x,y = 1, c z = −2, e the unit of charge, and ǫ 0 the electric constant. For a simulation of the Haldane-Shastry model, it turns out to be convenient to generate all interactions using transverse mode. Assuming isotropy in the transverse directions, we will be able to suppress the index α in the following.
It is convenient to transform the h α into the interaction picture of H 0 : h α → e iH0t/h h α e iH0t/h . This amounts for replacingâ m → a m e iωmt and σ x = 1 2 (e iω0t σ + ) + H.c. Finally, we make a rotating-wave approximation, that is, we neglet all fast oscillating terms in the Hamiltonian. To do this, we have to make some assumptions about the energies involved in the Hamiltonian: We tune the frequencies ω x and ω y , that is, the frequencies of the light field in h x and h y of Eq. (A2), close to ω 0 + ω m , that is, the frequency of a spin flip and the creation of a phononic mode. With this choice we obtain: Here, we have set the phase of the k x -laser field, φ x , set to zero. For the coupling in y direction, we choose this phase to be φ y = π, which effectively replaces σ x by the σ y matrix. Note that in the rotating-wave approximation, we have also neglected terms oscillating with ω α − ω 0 + ω m , which is justified if we tune ω α − ω 0 towards the upper edge of the phonon spectrum. Since no spin flip should be associated to the h z coupling, the corresponding frequency ω z must be tuned close to ω m . In the rotating-wave rotation, the corresponding Hamiltonian reads Unless ω mα = ω nβ , these functions only yield oscillatory terms, and can then be neglected. In that case, the time evolution of the system is equivalent to the one of a spin system with Hamiltonian