Complete devil's staircase and crystal--superfluid transitions in a dipolar XXZ spin chain: A trapped ion quantum simulation

Systems with long-range interactions show a variety of intriguing properties: they typically accommodate many meta-stable states, they can give rise to spontaneous formation of supersolids, and they can lead to counterintuitive thermodynamic behavior. However, the increased complexity that comes with long-range interactions strongly hinders theoretical studies. This makes a quantum simulator for long-range models highly desirable. Here, we show that a chain of trapped ions can be used to quantum simulate a one-dimensional model of hard-core bosons with dipolar off-site interaction and tunneling, equivalent to a dipolar XXZ spin-1/2 chain. We explore the rich phase diagram of this model in detail, employing perturbative mean-field theory, exact diagonalization, and quasiexact numerical techniques (density-matrix renormalization group and infinite time evolving block decimation). We find that the complete devil's staircase -- an infinite sequence of crystal states existing at vanishing tunneling -- spreads to a succession of lobes similar to the Mott-lobes found in Bose--Hubbard models. Investigating the melting of these crystal states at increased tunneling, we do not find (contrary to similar two-dimensional models) clear indications of supersolid behavior in the region around the melting transition. However, we find that inside the insulating lobes there are quasi-long range (algebraic) correlations, opposed to models with nearest-neighbor tunneling which show exponential decay of correlations.


I. INTRODUCTION
Dipolar interactions in a lattice system can lead to exciting new physics not present in systems with shortrange interactions [1]. For example, they can give rise to new ground states, like crystal states with fractional filling factors or supersolids. For instance, on the square lattice under hole-doping from half-filling, dipolar interactions can stabilize a supersolid phase [2], in contrast to interactions of limited range [3][4][5]. Further, dipolar interactions can lead to the appearance of a large number of metastable states [6]. On the one hand, these can complicate finding the ground state in numerical studies, while on the other hand they might be useful for applications in quantum information and as quantum memories [7][8][9]. Systems with strong long-range interactions (where the integral over the interactions does not converge) can even exhibit counterintuitive thermodynamic effects like super-extensive behavior, breaking of ergodicity, or a non-concave entropy [10]. These effects, however, are not the subject of the present work, since we consider weak long-range interactions.
Since long-range interactions can make theoretical calculations more difficult, it would be highly desirable to study all this extremely rich physics with an analog quantum simulator -specially for dynamics or two or more * Electronic address: Philipp.Hauke@icfo.es dimensions. In one dimension (which is the subject of this work), standard analytical and numerical techniques typically suffice to give a good picture of the phase diagram, but this is an important limit to gain intuition. Its experimental implementation allows to judge the reliability and performance of analogue quantum simulators. Natural candidate systems for this type of quantum simulation are trapped ions, where the high degree of control over state preparation, evolution, and readout allows access to a wide range of models. In 2004, trapped ions were first proposed for quantum simulations of certain classes of spin models by Porras and Cirac [11]. In a different context, similar spin models were derived earlier by Mintert and Wunderlich [12], who studied the possibility of individual ion addressing using inhomogeneous magnetic fields. Recent proof-of-principle experiments showed that, indeed, quantum spin systems can be simulated faithfully by systems of cold trapped ions [13,14]. For a review of trapped ion quantum simulations see [15].
An advantage of trapped ions is that dipolar interactions (which are the most common type of long-range interactions) are naturally achieved without additional experimental effort [11], since the interactions are mediated by phonons. This leads also to the peculiarity that, contrary to other possible setups for the quantum simulation of spin models, not only two-body (density-density) interactions, but also tunneling terms can be made long ranged.
In this work, we consider an experimentally feasible setup which allows to quantum simulate a chain of XXZ arXiv:1008.2945v2 [cond-mat.quant-gas] 3 Dec 2010 spins with dipolar interactions in a magnetic field. The system, which, as we show, supports a rich ground-state phase diagram, is described by the Hamiltonian where µ is the chemical potential (equivalent to an external magnetic field), and the S α i are spin-operators at site i.
The problem can be mapped to a system of hard-core bosons using the standard Holstein-Primakoff transfor- i (a i ) creates (destroys) a boson at site i and where n i denotes the corresponding number operator. In this picture, an up-spin corresponds to an occupied site and a down-spin to an empty one, the ZZ terms translate to dipolar off-site density-density interactions, and the XX and YY interactions become long-range tunneling terms. Since typically tunneling terms are short-ranged, the latter has to be seen as an important peculiarity of our model. In the following, we will use both the spin and boson pictures interchangeably, because some aspects are better described in terms of spins, while others are more familiar in the language of hard-core bosons.
By varying the angle θ, we can explore all ranges of relative strength of the interactions, e.g. the XX chain in a transverse field (θ = π/2), or the Ising model (θ = 0) -but both with long-range interaction. Moreover, the model can be tuned from negative to positive XY interaction, with the latter leading to a deformation of the phase diagram due to frustration effects.
At θ = 0, where Hamiltonian (1) behaves essentially classical, this model is known to show a succession of insulating states with different rational filling factors [16]. This succession, called 'devil's staircase', covers the entire range of the chemical potential, with (for an infinite system) each rational filling factor occuring on a finite range of chemical potential. Since the Hamiltonian (1) shows a symmetry of up-down spins (or in hard-core boson terminology, a particle-hole symmetry), we consider in this work only negative magnetization (filling factors lower than 1/2). When |θ| grows, the insulating states become destabilized and finally melt to a superfluid phase. We concentrate on characterizing the complete phase diagram of the model, and the nature of the correlations in each phase. We find that the long-range nature of the interactions changes qualitatively one important aspect of the physics of the problem: inside the insulating phases, the off-diagonal correlations decay algebraically instead of exponentially (as is typical in models with short-range tunneling) [17]. We call this phase a quasi-supersolid, since this is the closest analogue to a supersolid that can exist in one dimension: A supersolid is characterized by the coexistence of diagonal and off-diagonal long-range order (LRO). In 1D, however, off-diagonal LRO cannot occur [18], and the 'best' that can exist is the coexistence of diagonal LRO and off-diagonal quasi-LRO with algebraically decaying correlations. This quasi-supersolid is exceptional, since normally systems with dipolar interactions in 1D can be described by Luttinger theory [19,20]. In that case, diagonal correlations decay with the Luttinger parameter K, and off-diagonal correlations with 1/K. This means that if one of them decays slowly, the other one decays very fast. In contrast, in our model with long-range tunneling both diagonal and off-diagonal correlations decay slowly, and Luttinger theory is not a valid description.
Recently, a similar model with short-range tunneling has been investigated by a strong-coupling expansion in the context of ultracold dipolar atomic gases [21]. Placed in a one-dimensional optical lattice, in the limit of strong on-site interaction, such a system becomes similar to the one described by Hamiltonian (1), but with nearestneighbor (NN) instead of long-ranged tunneling. Hence, the quasi-supersolid phase cannot be observed in this system. Throughout the paper, we draw comparisons between our model and the model with NN tunneling and dipolar interaction, as well as with a model with both NN tunneling and interaction (the NN-XXZ model).
We organize this paper as follows: First, we briefly explain in Sec. II how the desired model Hamiltonian H can be implemented in a system of trapped ions. Then, in Sec. III, we outline our expectations by discussing previous related results. The following sections are dedicated to a thorough investigation of the ground state phase diagram of the model. Our most accurate analysis of the problem comes from numerical density matrix renormalization group (DMRG), Sec. V. However, it is good to first form intuition via analytical approaches, even if approximate. We will present mean field and perturbative results in Sec. IV to find the upper borders of the crystal lobes. Some limitations of DMRG can be overcome with other numerical techniques. For instance, we can study infinite systems (as opposed to finite sized) with the infinite time evolving block decimation (iTEBD) algorithm (section Sec. VI). Finally, experimentally relevant small systems can be thoroughly studied with exact diagonalization (Sec. VII). In Sec. VIII, we offer some final remarks.

II. EXPERIMENTAL IMPLEMENTATION
In this section, we shortly explain how the desired Hamiltonian, Eq. (1), can be simulated in an ion trap experiment. The discussion is a slight generalization of Ref. [11], to which we refer the reader for details. By applying a standing laser wave to a chain of trapped ions, one can arrive (after a canonical transformation) at an effective spin model where two internal hyper-fine states act as a pseudo-spin. In this scenario, the Coulomb interaction between two ions i and j transmits an effective spin-spin interaction where F α is the field strength of the laser propagating along direction α = x, y, z, m is the ion mass, and the elasticity matrix is Here c x,y = 1, c z = −2, and ω α is the vibration frequency of an individual ion. In the limit where the energy scale given by the Coulomb interaction between neighboring ions becomes much smaller than the energy scale of the vibration of the individual ions ω α , i.e.
1, the decay of the spin-spin interactions obeys a dipolar power law. Here, we introduced the mean inter-ion distance d 0 , and assumed that the ion chain is homogeneous. In a realistic experimental situation, if an overall trapping potential is used, the inter-ion spacing becomes non-uniform. In the central region of a long ion chain, however, the inhomogeneity can be neglected. Another way of creating an equidistant chain would consist of placing the ions in individual microtraps [22]. For simplicity, we will focus on the homogeneous situation.
Using three pairs of laser waves, one can thus induce the desired dipolar spin interactions in all the three directions of spin space. Choosing the detuning of the laser beams appropriately, one can tune the spin-spin coupling to negative or to positive, as desired. The missing chemical potential term can be achieved rather easily via an rf-frequency field driving the transition between the two hyper-fine states which consitute the effective spin 1/2. This completes all terms in Hamiltonian (1).

III. EXPECTED BEHAVIOR OF THE MODEL
Before proceeding to the detailed numerical analysis of the ground state phase diagram, let us briefly sketch the expected behavior of the present model, starting from previous related results.
For θ = 0, Hamiltonian (1) describes a 1D system of localized hard-core bosons with repulsive dipolar interaction, a classical model that can be solved analytically [16]. Due to the long-range nature of the interactions, for any given filling factor, the particles arrange in a periodic crystal pattern (a generalized Wigner lattice). For a given filling factor, these periodic patterns can be constructed by maximizing the mutual distance between the particles. Every rational filling factor q = m n occupies a finite extent in µ, thus giving rise to a plateau of fixed particle density, and these plateaux cover the entire range of µ. Plotting the filling factor against the chemical potential µ yields a self-similar structure -a complete devil's staircase, similar to Ref. [16]. The stair's steps have a width ∆µ m n = 2 3nπ 2 csc π n 2 − nπ 2 − 3π 3 cot π n csc π n 2 3n 3 .
(4) This means that fillings with a large denominatorequivalent to a large crystal period -become the ground state only in a very small range of chemical potential.
A finite tunneling, i.e. (θ mod π) = 0, introduces quantum mechanics into the system. The tunneling allows particles to gain kinetic energy, which destabilizes the crystal until it melts at some critical tunneling strength. This gives rise to a superfluid (SF) phase where the particles are delocalized over the chain.
There are two possible scenarios for a transition from the crystal structure to the molten phase: a direct crystal-SF transition, or an intervening supersolid phase. The latter is an exotic quantum phase showing crystal and SF behavior at the same time, which has been predicted for similar -but two-dimensional -hard-core boson systems [5,[23][24][25][26][27]. Until recently, there was only one claim of an experimental realization of a supersolid [28,29] in 4 He, which is still disputed [30,31]. Therefore, it would be interesting to find other systems which show supersolid behavior and are cleaner to interpret. One such experiment has been carried out very recently [32], where an atom cloud in a cavity breaks spatial symmetry while at the same time showing off-diagonal long-range order.
In the limit of θ = −π/2, Hamiltonian (1) describes a ferromagnetic XY model. Here, the most important effects of the long-range nature of the tunneling can be captured in a renormalized nearest-neighbor (NN) interaction: first, all interactions work towards the same ordering, and, second, in one dimension the integral over dipolar interactions converges. For θ = +π/2 the system is a frustrated XY antiferromagnet and the behavior is less obvious. It turns out, however, that a similar renormalization to a NN interaction captures the main physics. Hence, analogies to a NN-XY model can help understanding the behavior of the system near θ = ±π/2.
In the following, we study how the ground state phase diagram of Hamiltonian (1) develops as a function of θ. We start with two approximative methods to obtain qualitative insights.

IV. MEAN-FIELD APPROXIMATIONS
A. Perturbative mean-field theory A first rough understanding of the crystal-SF transition can be obtained by a perturbative mean-field theory (PMFT), valid for small tunneling. While not being very accurate in 1D, such a mean-field treatment generally becomes better for longer-ranged interactions.
Following the standard recipe of this Ansatz, we decompose the Hamiltonian into two parts, H = H 0 + H 1 . We choose i.e. we consider the insulating states, the exact ground states at the devil's staircase (θ = 0), as the unperturbed states. Afterwards, we introduce the perturbation The aim is now to calculate the expectation value of the SF order parameter where S − i is the spin lowering operator at site i and ρ is the density matrix of the system. At the tunneling strength where ϕ i becomes finite, the assumption that the ground state is localized is no longer valid. This gives an upper border for the crystal lobes.
For the evaluation of the expectation value in Eq. (7), we need the density matrix ρ, which is given by with Z the partition function and β = 1/ (k B T ) the inverse temperature. We are only interested in ground state properties, equivalent to the limit β → ∞. In this limit, Z → e −βE0 , where E 0 is the ground state energy. For small tunneling, i.e. small H 1 , we can approximate ρ using the Dyson expansion In mean-field approximation, we have (10) Inserting Eqs. (8) to (10) into (7), we obtain the following self-consistent equations for the SF order parameter where we introduced the abbreviation It can be checked that in the trace only one-hole excited states and the ground state itself have to be taken into account. For small J tan θ, Eq. (11) has only the trivial solution ϕ i = 0. At some critical J c (µ), however, SF order may develop, i.e. ϕ i = 0. This is equivalent to an instability under adding or removing a single particle. The point where this happens gives an upper border for the crystal phase. The main disadvantage of the current method is that there are potentially more complicated excitations, e.g. the addition of a particle plus a relocation of the neighboring particles. Such a deformation of the crystal could decrease the potential energy. In the current simple Ansatz, however, this type of excitations is not captured. For our calculations, we assumed an infinite system with the restriction that the states have a periodicity of 12 sites. Since the crystal phases with a large periodicity become very small in their extent in µ [by virtue of Eq. (4)], this is a reasonable restriction which still captures the most prominent features of the phase diagram.
The insulating-lobe structure that we obtain is shown in Fig. 1. The thick line follows the breakdown of the ground state. A dotted line marks the values of µ where the ground state magnetization changes, i.e. where a different crystal structure becomes lower in energy. However, even when a given crystal structure ceases to be the ground state, it can still be a metastable state with respect to one-particle or one-hole excitations. In Fig. 1, thin lines mark the corresponding regions where the most important crystal states are metastable. For clarity, states that constitute the ground state over only a very small region of the phase diagram have been excluded. Figure 1 shows the rich structure of metastable crystal states. Recently -in the context of ultra-cold dipolar neutral atoms -it has been proposed to use such metastable states as quantum memories [8,9]. Another feature of the crystal structure is that it is more stable for frustrated antiferromagnetic XY-interaction (θ>0) than for ferromagnetic XY-interaction (θ<0).
We also computed the phase diagram in a Gutzwiller mean-field theory (not shown). The qualitative behavior is similar, but the lobes are somewhat smaller. The reason is that PMFT considers destabilization under singleparticle or single-hole excitations, while Gutzwiller meanfield theory captures better more complicated excitations.

B. Wigner-crystal melting at low filling
For very low magnetizations, we can draw an analogy to the melting of Wigner crystals [33,34]. As described in Sec. III, at θ = 0, the hard-core bosons are perfectly localized, with an inter-particle distance given by the filling fraction. At finite tunneling, θ = 0, the particles spread, but at small tunneling and low filling it is reasonable to assume that they remain spatially well separated. Under this assumption, we can approximate the   total wave-function by a product of Gaussians representing the individual particles [43]. Using this Ansatz, one can obtain a self-consistent equation for the distribution of any particle over the lattice.
The spread of the wave packets increases with θ until the initial assumption that the particles are well separated breaks down, and the crystal has to be considered as molten. A measure for this is the normalized variance L ≡ ∆x d0 = 1 d0 x 2 − x 2 , called the Lindemann parameter [35]. As shown in Fig. 2, for large inter-particle spacing d 0 , the system behaves similar to the continuum limit, where it is known [33] that there is a melting transition at J tan θ c = const/d 0 . We find that the contour lines of the Lindemann parameter from our crude approximation follow this behavior well for large d 0 . At small inter-particle distances there are deviations, but here the For dipolar hopping, the decay is algebraic everywhere, even within the lobes (A and C), where it follows the dipolar exponent α = 3, while for NN hopping the decay is exponential. In the SF (B and D), the algebraic decay is faster than for NN hopping on the antiferromagnetic side and slower on the ferromagnetic side. Panel E shows a zoom on the dashed section together with the finite-size scaling result of the crystal lobe at 1/3 filling, as explained in the text. assumptions made in our approximation are not valid anyway. The physical explanation behind the observed behavior is that for large inter-particle distances the repulsive interaction becomes very weak. Hence, for larger d 0 a smaller gain in kinetic energy suffices to destabilize the crystal.

V. DENSITY-MATRIX RENORMALIZATION GROUP (DMRG)
While mean-field theories can provide some physical understanding of the system, in 1D they are not very precise. Therefore, we turn now to a thorough numerical analysis of the phase diagram by the quasi-exact densitymatrix renormalization group (DMRG) method [36][37][38]. We consider chains with open boundary conditions of up to N = 102 sites, and with a bond dimension of D = 128. The interaction range was always cut-off at half the length of the chain. Figure 3 presents the mean polarization Z /N ≡ 1 2 N j=1 σ z j /N , where N is the number of sites and σ z j the Pauli-z matrix of spin j. We find that for large enough field µ, or strong ferromagnetic ZZ-interaction (corresponding to θ not too far from ±π), the system is in a fully polarized state, quantitatively similar to what we saw in the mean-field calculations. However, as expected, the general precision of the mean-field calculations in this one-dimensional system is low: the DMRG results indicate that the θ-range of the crystal lobes is up to an order of magnitude smaller than predicted in meanfield. In fact, in the global view given in the main panel of Fig. 3, only the plateau for 1/2 filling is discernible. Panel E of Fig. 3 shows the small region of the phase diagram where the 1/3-filling crystal lobe is located. Due to the long-range nature of the hopping terms, the polarization is asymmetric with respect to ±θ, in contrast to [21].
We find that open boundary conditions play an important role for the finite systems used here. For example, in the center of the main panel of Fig. 3, there appears a broad plateau with (1/2 − 1/N ) filling (i.e. with one spin flipped away from half filling). In the thermodynamic limit, this phase is not as stable and merges with the 1/2filling plateau to a single lobe. This peculiarity is highly relevant for the correct interpretation of future experiments, which might be carried out with open boundary conditions for technical convenience. Moreover, boundary effects play an important role in the determination of the decay of correlations, which we now turn to describe.
The log-log plots in panels A to D show the decay of the in-plane correlations G (j) ≡ σ + N/2 σ − N/2+j at some selected points in the phase diagram (filled circles), and compare the result to the same system but with nearestneighbor (NN) hopping (solid line). In the superfluid phase (panels B and D), the decay is algebraic for both NN and dipolar hopping. Panels A and C lie inside the lobes corresponding to 1/2 and 1/3 filling, respectively, which is visible through the oscillation of the correlation function superposed with an overall decay. We see that, for dipolar hopping, the decay is algebraic everywhere, even within the lobes, where for NN hopping the correlations decay exponentially. In the lobes, the decay for dipolar hopping follows the interactions with exponent α = 3 [17]. Notice that a clear fit for this exponent can only be obtained for rather large chains with N > 60 spins. Even more, although the strength of the XY-interactions decays very rapidly with distance, the power-law decay of correlations inside the lobes cannot be observed at all if the interaction is truncated.
Since inside the lobes we have a crystal structure, the long-range decay of the transverse correlations is an anomaly when compared to other models with insulating phases. For this similarity to a supersolid, we call this phase a quasi-supersolid phase. In a true supersolid, a crystal structure coexists with a long-range order in the transverse correlations. However, in 1D off-diagonal LRO cannot be spontaneously broken [18]. Hence, the closest analogue to a supersolid that can exist in 1D is a state as the present one -with diagonal LRO and algebraically decaying off-diagonal correlations.
A change in the behavior of the correlations marks the transition from a crystal state to the SF. For a quantitative evaluation, at each point in the phase diagram we fit an algebraic trial function G a (j) = c 0 j −α . The value of α for a chain of N = 60 spins is shown in Fig. 4. The cuts of Fig. 4 at µ = 0 (1/2 filling, A) and θ = π/2 (XY model, B) show the wide range that the exponent of the correlations takes in this system. They also demonstrate the influence of frustration: the correlations decay faster at the antiferromagnetic side (θ > 0) than at the ferromagnetic side (θ < 0). The exponents are similar to the NN XXZ model (black line in A), but cover a broader range. In particular, they are not independent of the filling for θ = π/2 (panel B), where the NN XXZ model has α = −0.5. The results for a system with dipolar ZZ-interaction but NN XY-interaction are shown in red (dashed line in A, filled circles in B). The main qualitative difference is that inside the insulating lobes the decay is no longer a power law, but follows an exponential. In a finite system, finding the phase transition from the crystal state to the SF is not simple, because the finite number of spins prevents the system from assuming arbitrary polarizations. This results in a division of the phase diagram into different stripes with fixed integer number of up spins, which makes the crystal lobes difficult to discern. In infinite systems, however, one expects a steplike behavior of the polarization only in crystal phases, while in the SF the polarization should change smoothly with µ between −1/2 and 1/2. This observation makes it possible to extrapolate the border between the crystal phases and the SF by the following finite size scaling: At fixed θ and for a given polarization, compute for several chain lengths N the upper and lower limits of the polarization plateau, µ a (N ) and µ b (N ). Then, from a finite size scaling of the results, one can extract µ a,b (N = ∞). In the SF, where the polarization should be a continuum, µ a (∞) and µ b (∞) should be equal. However, in the crystal lobe there will be a finite distance between µ a (∞) and µ b (∞): the width of the lobe for that value of θ. A similar procedure is known to work well for the estimation of Mott-lobes from finite systems in the Bose-Hubbard model [39]. The main difference is that in our model there is in principle a lobe for any rational filling factor, instead of only for integer filling factors. Panel E of Fig. 3 shows the result of this approach for the 1/3filling lobe, using chain lengths of up to 102 spins. The cusp structure is typical for one dimensional systems.

VI. INFINITE TIME EVOLVING BLOCK DECIMATION (ITEBD)
To complement the results of the previous section, we study the phase diagram with the iTEBD algorithm [40], in which infinite chains can be directly addressed without the need of finite size extrapolations. Treating long-range interaction complicates the original formulation of the iTEBD algorithm and turns out to lead to convergence problems. Instead, we implement a variant of the original algorithm in which the translationally invariant character of the problem is kept in the state and the interactions, thanks to the use of Matrix Product Operators [41]. In this way, we can include interaction terms ranging longer than nearest neighbors in a much simpler way.
By comparing the phase diagrams for different range of XY interactions, we may study the effect of the longrange hopping in the thermodynamic limit. The method still requires a truncation of the interacion range to some finite order, so that it will not be possible to reproduce the power-law decay of correlations within the crystal lobes.
We truncate the ZZ interactions above next-to-nearestneighbor (NNN) interactions and compare the cases of NN and up to NNN XX and YY terms. We observe that a small bond dimension, D = 10, provides already a good approximation to the overall phase diagram. To analyze the isolating lobe at 1/3 filling we increase the bond dimension to D = 20. Figure 5 shows the results for the magnetization per particle in the first case, in which XX and YY interactions range only to the NN. As in [21] the phase diagram is symmetric. This is also seen in the close up on the filling 1/3 crystal lobe (right panel of Fig. 5). Figure 6 shows the equivalent phase diagram when NNN terms are included in the XY plane. Including one more term in the XY interactions already causes a clear asymmetry of the phase diagram with respect to the change θ → −θ. This is clearly visible in the superfluid phases (left panel of the figure), and the zoom on the isolating lobe at polarization −1/6 shows some asymmetry as well.
The size of the lobe is larger than the result from DMRG after finite size extrapolation. We are, however, only including the first term further than the NN, and we expect that longer range terms correct the exact shape of the lobe, as they will give rise to lobes corresponding to other filling factors.

VII. EXACT DIAGONALIZATION
In this section, we supplement the previous results by exact diagonalization (ED) of chains of 18 spins with periodic boundary conditions, using the Lanczos method. Although ED only allows investigation of very small systems, from these one can often infer surprisingly much about the behavior of the model at larger sizes. More- over, an experimental implementation with a chain of trapped ions will have to start with short chain lengths. For such a case ED provides a valuable validation of the trapped-ion quantum simulator. Figure 7 shows the fidelity susceptibility, for 1/3 filling. At a second order quantum phase transition, χ F diverges in the thermodynamic limit. The peaks in a finite system are precursors thereof. However, in the present system the divergence turns out to be relatively weak, which makes it difficult to extrapolate quantitative values for the critical point of the transition. The qualitative value, however, is consistent with the DMRG and iTEBD analyses. The second derivative of the energy has been proposed [42] as a substitute of the fidelity susceptibility to detect quantum phase transitions. In the present system, however, the energy is very smooth and does not give a good indication of the phase transition. Instead, we compare in Fig. 7 χ F to the second derivative of the total in-plane magnetization, Here, the prime indicates that the sum excludes self-correlations. We take the absolute value of the correlations in order to take into account the different staggered configurations occurring for different filling factors. We see that ∂ 2 M x /∂θ 2 peaks at slightly larger absolute values of θ than χ F , but the rough locations are consistent with the peaks of χ F [44].
A comparison to the results from DMRG shows a very similar qualitative behavior. This means that already the very small chains considered here show strong precursors of the physics relevant for large chain lengths.

VIII. CONCLUSION
To summarize, we have studied a one dimensional system of hard-core bosons where particles interact and can tunnel at long distances with an algebraically decaying strength. We have discussed how the system can be mapped onto a spin Hamiltonian that can be simulated experimentally using trapped ions. Unlike other atomic quantum simulators (such as dipolar ultracold atoms), trapped ions appear to be more flexible in the manipulation of some parameters -e.g. it is easy to change the interactions and hoppings from negative to positive simply by detuning the lasers.
The system we have studied has a rich phase diagram, with many insulating phases filled with (possibly longlived) meta-stable states and quasi-long range order in all correlations -two prominent features that are induced solely by the long-range nature of the interactions.
The coexistence of diagonal LRO and off-diagonal quasi-LRO in the insulating lobes can be seen as a quasisupersolid, the closest analogue to a supersolid that can be found in 1D. It would be very interesting to study how these phases behave in 2D, where the possibility of having off-diagonal LRO could allow the quasi-supersolids to become true supersolids.
Since calculations prove to be quite complicated for long-range interactions, one can consider the present model as a testbed for a comparison of different numerical tools. We find that the precision of the standard perturbative mean-field theory is insufficient, because it does not capture excitations more complicated than single particle excitations. Exact diagonalization, DMRG, and iTEBD, however, all deliver a consistent picture. Among these, we find that DMRG is the method of choice for the model studied. However, it is restricted to finite systems. While iTEBD is less accurate and restricted to shorter range of interactions, it allows to investigate the thermodynamic limit directly. Exact diagonalization, finally, while being relevant only for very small systems, can still serve for benchmarking other less accurate computational methods and near-future experiments.
In the future, it would be interesting to pursue investigations of other properties of the system that are affected by the long-range nature of the interactions, such as its response to excitations or the dynamics of its correlations.