Quantum spin models with long-range interactions and tunnelings: a quantum Monte Carlo study

Michał Maik; Philipp Hauke; Omjyoti Dutta; Jakub Zakrzewski; Maciej Lewenstein

doi:10.1088/1367-2630/14/11/113006

1. Introduction

Quantum simulators (first proposed by Feynman [1]), which are devices built to evolve according to a postulated quantum Hamiltonian and thus 'compute' its properties, are one of the hot ideas which may provide a breakthrough in many-body physics. While one must be aware of possible difficulties (see, e.g., [2]), impressive progress has been achieved in recent years in different systems employing cold atoms and molecules, nuclear magnetic resonance, superconducting qubits and ions. The latter are extremely well controlled, and it has already been demonstrated that, indeed, quantum spin systems may be simulated with cold-ion setups [3, 4]. Quantum spin systems on a lattice constitute one of the most relevant cases for quantum simulations as there are many instances where standard numerical techniques for computing their dynamics or even their static behavior fail, especially for two-dimensional (2D) or 3D systems.

Chronologically, the first proposition to use trapped ions for simulating lattice spin models came from Porras and Cirac [5], who derived the effective spin Hamiltonian for the system⁶,

$\begin{equation} H = J\sum\limits_{i,j}\frac{1}{|i-j|^{3}}[\cos{\theta}(S_{i}^{z}S_{j}^{z})+\sin{\theta}(S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y})]-\mu\sum\limits_{i} S_{i}^{z}, \end{equation} \tag{ 1 }$

where μ is the chemical potential (which, in this case, acts as an external magnetic field), J > 0 is the interaction strength and S^α_i are the spin operators at site i. All the long-ranged interactions fall off with a 1/r³ dipolar decay, which for the ions is due to the fact that they are both induced by the same mechanism, namely lattice vibrations mediated by the Coulomb force [7]. Hamiltonian (1) is a dipolar XXZ model where the ratio of hopping to dipolar repulsion can be scaled by θ [8, 9]. Another route to simulate quantum magnetism of dipolar systems, using the rotational structure of the ultracold polar molecules, has been proposed in [10, 11]. In this case, the XY (or tunneling) term is restricted to nearest neighbors (NN), and only the ZZ (or interaction) term is dipolar.

Precisely this long-range (LR) dipolar character makes this system interesting and challenging. Dipolar interactions introduce new physics into conventional short-range (SR) systems. For example, for soft-core bosons, the system with LR interactions and SR tunneling can—apart from Mott-insulating, superfluid and crystal phases—host a Haldane-insulating phase [12]. This is characterized by antiferromagnetic order between empty sites and sites with double occupancy, with an arbitrarily long string of sites with unit occupancy in between⁷. Another case when dipolar interactions play a crucial role is the appearance of the celebrated supersolid. The extended Bose–Hubbard model (i.e. with NN interactions) with soft-core on-site interactions shows stable supersolidity in one-dimensional (1D) [13] and square lattices [14], where in the 1D case a continuous transition to the supersolid phase occurs, in contrast to the first-order phase transition in 2D lattices. Long-ranged dipolar interactions, on the other hand, allow for the appearance of supersolidity even in the hard-core limit [15, 16]. In this limit, dipolar interactions give rise to a large number of metastable states [17, 18] (for a review see [19]). By tuning the direction of the dipoles, incompressible regions like devil's staircase structures were predicted in [20]. Despite being interesting, LR interactions make computer simulations of such systems quite difficult. On the other hand, since ions in optical lattices may be extremely well controlled, they form an ideal medium for a quantum simulator.

In our previous paper [8], we considered both the mean-field phase diagram for the system (1) as well as a 1D chain using different quasi-exact techniques, such as the density matrix renormalization group, and exact diagonalization for small systems. Here, we would like to concentrate on frustration effects on a 2D triangular lattice using a quantum Monte Carlo (QMC) approach.

For the system consisting of spin-half particles (such as, e.g., originating from two internal ion states in the Porras–Cirac [7] proposal), the spins can be mapped onto a system of hard-core bosons, using the Holstein–Primakoff transformation S^z_i → n_i − 1/2, S⁺_i → a^†_i and S⁻_i → a_i, where a^†_i (a_i) are the creation (annihilation) operators of hard-core bosons and n_i is the number operator at site i. These bosons obey the normal bosonic commutation relationships, [a,a^†] = 1, but are constrained to only single filling, a² = (a^†)² = 0. This is the limit when the on-site repulsion term U goes to infinity in the standard Bose–Hubbard model. In this representation, a spin-up particle is represented by a filled site and a spin-down particle by an empty site. For the sake of convenience, we will mostly use the language of hard-core bosons in this paper.

Before presenting the results for the 2D triangular lattice, let us briefly review the behavior of the system for a 1D chain of bosons as obtained in [8].

2. Review of the results in one dimension

In order to form some intuition about possible physical effects due to the LR interaction and tunneling, we now discuss briefly the ground-state phase diagram of the Hamiltonian in 1D. For hard-core bosons, which is the topic of this paper, the 1D version of the Hamiltonian has been thoroughly investigated in the past [8] (see also [21, 22] for the special cases θ = 0 and θ = π/2). For reference, we reproduce in figure 1 the phase diagrams for the system with dipolar interaction and NN tunneling, as well as the system with both the interaction and the tunneling terms dipolar.

**Figure 1.** Phase diagram of the 1D system with dipolar interactions and (a) NN tunneling and (b) dipolar tunneling. The phases are labeled according to density matrix renormalization results, while the actual data shown come from the infinite time evolving block decimation method (with interactions truncated at next-to-NN), both from [8]. Along the black line, there exists a devil's staircase of crystal states. At finite tunneling, for (a), these spread into conventional insulating states, while for (b), they form quasi-supersolids. The light blue lines near the center of the plots sketch the crystal lobes at 2/3-filling. Their cusp-like structure is typical of 1D systems. In 2D, they are expected to be rounded off, similar to Mott lobes of the Bose–Hubbard model [23]. For NN tunneling (a), the superfluid (SF) phases can be mapped into one another, while for NN dipolar tunneling (b), they are distinct on the ferromagnetic (FM, θ < 0) and the antiferromagnetic side (AFM, θ > 0). Note also that in (b), frustration leads to an asymmetry between θ < 0 and θ > 0.
Download figure:
Standard image

At zero tunneling, the ground states are periodic crystals where—to minimize the dipolar interaction energy—occupied sites are as far apart as possible for a given filling factor [21]. For finite 1D systems and very small tunneling, such a situation persists as exemplified in [24]. For infinite chains in 1D, every fractional filling factor n = p/q is a stable ground state for a portion of the μ parameter space. The extent in μ decreases with q, since at large distances the dipolar repulsion is weak and thus cannot efficiently stabilize crystals with a large period. This succession of crystal states is termed the devil's staircase. This name derives from its surprising mathematical properties, challenging naive intuitions about continuity and measure: since all rational fillings are present, it is a continuous function; moreover, its derivative vanishes almost everywhere (i.e. it is non-zero only on a set of measure zero)—and still it is not a constant, but covers a finite range.

At finite tunneling, the crystals spread into lobes similar to the Mott lobes of the Bose–Hubbard model. If the tunneling is only over NNs, these Mott lobes are not sensitive to the sign of the tunneling and form standard insulating states with diagonal long-range order (LRO) and off-diagonal SR order (figure 1(a)). For LR dipolar tunneling, the extent of the lobes is asymmetric under sign change θ → − θ: frustration effects stabilize the crystal states for θ < 0, while the ferromagnetic tunneling for θ > 0 destabilizes them (figure 1(b)). Moreover, the crystal states acquire off-diagonal correlations which follow the algebraic decay of the dipolar interactions [8]. This coexistence of diagonal and off-diagonal (quasi-)LRO turns the crystal states into (quasi-)supersolids. The occurrence of such phases for hard-core bosons in 1D is truly exceptional, since the systems where it appears consist typically of soft-core bosons [13, 14] or 2D lattices [15]. Furthermore, this 1D quasi-supersolid defies Luttinger-liquid theory, which typically describes 1D systems very well, even in the presence of dipolar interactions [25–27]. Where Luttinger-liquid theory applies, diagonal and off-diagonal correlations decay algebraically with exponents which are the inverse of one another. Therefore, if the diagonal correlations show LRO, the corresponding exponent is effectively 0, and the exponent for the off-diagonal correlations is infinite, describing an exponential decay. In our case, this exponent remains finite in the quasi-supersolid phase and the above relationship clearly does not hold.

At even stronger tunneling strengths, the crystal melts and the system is in a superfluid phase. The LR tunneling and interactions influence the correlation functions at large distances and therefore also modify the character of this phase [8, 22].

These results show that in this system the dipolar interactions considerably modify the quantum-mechanical phase diagram. In higher dimensions, we can expect the influence of LR interactions to be even stronger, which makes extending these studies to a 2D lattice highly relevant. For example, one can expect that—if quasi-supersolids appear already in 1D—the LR tunneling has a profound effect on the stability of 2D supersolids, which appear in triangular lattices at the transition between crystal and superfluid phases [13, 28–30]. Also, the frustration effects already observed in the 1D system should be much more pronounced in the triangular lattice, simply due to the increased number of interactions⁸.

Further, such an analysis is especially relevant at finite temperature. In fact, a recent work scanned the phase diagram of Hamiltonian (1) along the line μ = 0 in a square lattice [9]. There, the authors found that above the superfluid on the ferromagnetic side (i.e. θ < 0), the continuous U(1) symmetry of the off-diagonal correlations remains broken even at finite temperatures. Thus, the LR nature of the tunneling leads to a phase which defies the Mermin–Wagner theorem [31]⁹.

For these reasons, we can expect intriguing effects of LR tunneling for the 2D triangular lattice, to which we now turn our attention.

3. Dipolar hard-core bosons on the triangular lattice

To analyze the system of dipolar bosons on the 2D triangular lattice, we employ QMC simulations at a finite but low temperature. Specifically, we study rhombic lattices with periodic boundary conditions and N = L × L sites, with L = 6–12. We will, in particular, thoroughly investigate the Wigner crystal at 2/3-filling. The triangular lattice is frustrated even with only NN interactions. When any longer-ranged interactions are added, this only increases the frustration, which also making QMC calculations more difficult. Also, in frustrated systems [32] and systems with LR interactions [18], typically many metastable states appear, rendering finding the ground state a somewhat difficult task. To avoid complications with the sign problem, caused by negative probabilities in QMC codes, we will take only negative values of θ into consideration. Note that this sign problem appears only for the long-ranged XY interactions, while the ZZ (Ising-like) interactions are—despite frustration—sign-problem free.

In this study, we compare the Hamiltonian with both LR interactions and hoppings (LR–LR), with long-ranged dipolar interactions but NN hopping (LR–SR; relevant for polar molecules [33]), as well as with hopping and interactions truncated at NNs (SR–SR; this is the NN XXZ model, relevant to magnetic materials with planar anisotropy in their couplings)¹⁰. Each of these systems will display different crystal, superfluid and supersolid regions. Comparing these cases will give valuable insight into the influence of the LR terms.

All the calculations were performed using the worm algorithm of the open source ALPS (Algorithms and Libraries for Physics Simulations) project [34]. This algorithm, first created by Prokof'ev, works by sampling world lines in the path integral representation of the partition function in the grand canonical ensemble [35].

3.1. Vanishing tunneling

The first calculation, which creates the motivation for the rest of the paper, is to look at the case of vanishing tunneling and temperature for each system (LR–LR, LR–SR and SR–SR). Here, similar to the 1D devil's staircase, at vanishing temperature a series of insulating crystal states is expected to cover the entire range of μ/J. Since we are interested in finite temperature results, we set T = 0.1—which should still be low enough to reflect the characteristics of the ground-state phase diagram—and look for plateaus in the density. We distinguish short- and long-ranged ZZ interactions. From figure 2, left panel, we can see that the only plateau (besides the completely filled system) that appears is at ρ = 2/3 (corresponding either to 2/3 boson filling or in spin terms, a lattice with 2/3 of the spins oriented up and 1/3 oriented down) for both short- and long-ranged interactions. Scaling the system size from L = 6 to 12 causes no change for the short-ranged interactions, and a minimal change for the long-ranged interactions. The key difference is in the size and position of the short-ranged and long-ranged plateaus. For short-ranged interactions this plateau is larger and centered around μ/J ≃ 1.5, while the long-ranged interactions have a smaller plateau centered around μ/J ≃ 1.85. The finite width of these plateaus suggests that a 2/3-filling Wigner crystal persists also for some finite θ. The right panels of figure 2 show how the plateau shrinks with increasing temperature, for SR interactions (top right panel) as well as for LR interactions (bottom right panel). In the latter case, in fact, by T = 0.25 the plateau has completely disappeared. We can also note that at T = 0.05 there are signs of some of the other plateaus, most noticeably the 3/4-filling plateau. The rest of the paper will focus on the 2/3-filling crystal lobes and their properties.

**Figure 2.** The graph on the left displays θ = 0 and T = 0.1, where the density shows a single plateau for 2/3-filling. For SR interactions (solid blue line), different curves for L = 6,9,12 coincide, and for LR interactions (solid red: L = 6; dashed red: L = 9 and 12) the size dependence is small. The panels on the right are at fixed L = 12 and θ = 0 for different temperatures T = 0.05,0.15 and 0.25 (from dotted to dashed to solid). Both for SR (top right) and LR interactions (bottom right), the 2/3-filling plateau shrinks as T increases.
Download figure:
Standard image

3.2. Low-temperature phase diagram at finite tunneling

We now introduce a finite tunneling by choosing θ < 0 in our Hamiltonian, and study the properties around the 2/3-filling crystal. We calculate the density and the superfluid fraction. The superfluidity is measured using the winding numbers calculated from the movement of the worms in the QMC code. In order to get this value the system must have periodic boundary conditions so that the world lines can properly 'wind' around the system. The superfluid fraction is

$\begin{equation} \rho_{\rm s}=\frac{\langle W^{2} \rangle }{4 \beta}, \end{equation} \tag{ 2 }$

where W is the winding number fluctuation of the world lines and β is the inverse temperature.

Figure 3 shows the results for the boson density of an L = 6 triangular lattice at T = 0.1. For all types of interactions, we see that as θ increases in absolute value, the $\rho = \frac {2}{3}$ plateau shrinks, because for larger |θ| the ratio of hopping to dipolar interactions increases. This introduces more kinetic energy and the crystal melts into a superfluid. It can be seen that for the SR–SR system the boson density lobe extends to θ ≃ − 0.36, while for the LR–SR interactions it ends at θ ≃ − 0.3, and finally for LR–LR interactions the lobe is still smaller, only ranging up to θ ≃ − 0.2. The behavior is explained by the fact that the increased number of interactions causes a quicker melting of the lobe.

A second major observation is the shift in the position of the lobes. While the short-ranged lobe exists approximately for 0.3 < μ/J < 2.7, the long-ranged lobes lie generally in the interval 1.0 < μ/J < 2.8. For a system like this at T = 0, one expects the $\rho = \frac {2}{3}$ lobe to exist on the range 0 < μ/J < 3 with a mirrored image of the $\rho = \frac {1}{3}$ lobe at −3 < μ/J < 0. The existence of the two identical lobes is explained by particle–hole symmetry [28–30]. The lobes are separated with a kind of mixed solid in between (with the coexistence of 1/3- and 2/3-filling regions). The existence of this region may be caused by several phenomena. It could be either an effect of the finite temperature as in [36] or due to the existence of many metastable states caused by a devil's staircase-like behavior, similar to what was observed in [20].

Again referring to the T = 0 phase diagram, we expect that there is a region of supersolidity that extends in between the lobes and goes all the way to their base at θ = μ = 0 [28]. In our system this region should exist near the tip of the lobe but not extend all the way to the base due to the finite temperature and the resulting mixed solid. Looking at figure 3, it is obvious indeed that, if a supersolid region exists, it can only be near the tip of the lobe because the superfluidity is zero a significant way up the lobe. Judging by the increased separation of the long-ranged lobes, we can assume that the supersolid region for these systems should increase in size to fill the region in between. To search for the supersolid phase, we now compare the superfluidity with the static structure factor.

The structure factor is defined as the Fourier transform of the density–density correlations,

$\begin{equation} S(\mathbf{Q})=\left\langle\left\vert \sum\limits_{i=1}^{N}n_{i}\,\mathrm{e}^{\mathrm{i} \mathbf{Qr}_{i}} \right\vert^{2}\right\rangle\Bigg/N^{2}. \end{equation} \tag{ 3 }$

Here, we focus on the wave vector Q = (4π/3,0), which corresponds to the $\sqrt {3} \times \sqrt {3}$ -order parameter that is associated with 1/3- and 2/3-filling crystals on the triangular lattice. For the case of the 2/3-filling lobe that we are interested in, it will show plateaus over the same range of μ as the density, but additionally gives insight into the arrangement of the bosons on the lattice. This makes it a useful quantity in searching for supersolid regions. In fact, a supersolid exists when both the structure factor and the superfluid fraction have non-zero values. The physical mechanism behind the supersolid phenomenon is based on the appearance of extra holes (particles). The underlying crystal structure has $\sqrt {3} \times \sqrt {3}$ order on a triangular sublattice of the physical lattice. The extra holes (particles) are free to move around in the rest of the lattice as superfluid objects. In this way, the system retains a crystal structure, while it acquires at the same time the LR coherence of a superfluid. Due to the hole (particle) doping, it forms in sections away from commensurate filling, in this case in between the 1/3- and 2/3-filling lobes.

Taking 'slices' out of the crystal lobes we now check where there is a supersolid region and where the system directly transitions from crystal to superfluid. Also we study the nature of these transitions to see if they are of first or second order.

The most logical region to look for supersolids is directly in between the 1/3- and 2/3-filling lobes, at μ/J = 0. Figure 4 shows the behavior of the SR–SR, the LR–SR and the LR–LR system at this cut for several system sizes. The structure factor and superfluidity reveal, for all the systems, three different regions. In each case, the system starts at θ = 0 in a solid phase where the superfluid fraction is zero but the structure factor is finite. It transitions smoothly into a supersolid region where both superfluidity and structure factor are non-zero. Finally, the structure factor smoothly drops away and leaves just a non-zero superfluid fraction, making the final phase a superfluid. In each system, the size of the supersolid region is different. In the SR–SR case, the supersolid region begins to appear at θ ≃ − 0.15 for an L = 6 lattice. As the size grows to L = 12, the region has shifted to θ ≃ − 0.19 with the superfluid curves becoming sharper. The increased system size also reduces the value of the structure factor a little. From [29], we know that at even larger sizes (but T = 0) the supersolid will continue to exist in this type of system. For the LR–SR system, the supersolid region appears at a similar point and also shifts with a system-size increase. The structure factor, on the other hand, has a significant decrease for larger system sizes. It is difficult to tell if at greater sizes the existence of the supersolid will persist. The final graph shows the LR–LR system. In this system, superfluidity appears even before μ/J reaches −0.1. In this system, the superfluid fraction is much greater than in the previous two because of the long-ranged tunneling. This means that at small system sizes the supersolid region is much more prominent relative to the crystal lobe. The structure factor diminishes with system size almost exactly as in the LR–SR case except that the transition is at a different value of θ, and near μ/J = 0 it drops to slightly lower values. Due to this strong decrease, for any situation with LR interactions we cannot clearly state whether the supersolid region survives at larger system sizes.

Next we take vertical cuts at a value of θ = −0.15, since this is a reasonable place for a supersolid to exist for the LR–LR system (≃80% of the tip of the lobe). We compare all the systems at this cut, and study the behavior of the superfluid fraction and the structure factor, plotted in figure 5. For SR–SR interactions, no supersolid region appears. On one side of the lobe there is a sharp phase transition directly from the crystal to the superfluid phase, while on the other side there is a slower change from one solid form to another (ρ = 1/3 → 2/3). The finite-size scaling in the figure shows that as the size increases the transitions of the structure factor become even sharper, although they stay continuous due to the finite temperature. The values of the superfluid fraction decrease as the system sizes increase and essentially disappear at L = 12. In the LR–SR system, a hint of the supersolid phase begins to appear on either side of the lobe. It is a bit more evident on the side where μ/J is small (as is to be expected from references such as [28]), but it also arises on the opposite side. This is contrary to a system of only short-ranged interactions where this supersolid region appears only on one side of the lobe and not both. At larger sizes, also in the LR–SR system the transitions become sharper and the superfluidity gets smaller. The final and most interesting cut is taken out of the LR–LR lobe. In this system, we see a smooth transition from crystal to supersolid at μ/J ≃ 2.4 and at μ/J ≃ 1.4 for L = 6. For larger systems the transition at μ/J ≃ 2.4 occurs at the same spot but becomes sharper, making the supersolid region disappear. At μ/J ≃ 1.4 the transition shifts to a higher value of μ/J, making the 2/3-filling plateau smaller. It also becomes less smooth but the supersolid region remains longer than for the SR–SR case. In the LR–LR system, the superfluid fraction for larger systems has the opposite effect than for the previous cases, it becomes higher instead of lower.

**Figure 5.** Cuts at θ = −0.15 for SR–SR (top left), LR–SR (bottom left), LR–LR (right) for a L = 6,9 and 12 lattice (the lines become thicker and darker with increasing system size). Solid blue: structure factor. Dashed black: superfluid fraction. At this value of θ, for SR tunneling, the superfluid fraction disappears rapidly with increasing system size, while for LR tunneling it even increases. At μ/J ≈ 1.75, possibly a supersolid may survive in large lattices. Again in the LR–LR system there is a reversal of finite-size effects. The superfluid fraction for larger systems becomes higher instead of lower.
Download figure:
Standard image

A perhaps fairer comparison is to look at a cut through a region where we are sure the supersolid exists for all three systems. Therefore, in figure 6 we look at two more cuts that are now taken closer to the tips of the SR–SR and LR–SR lobes. As with the LR–LR system (the last panel of figure 5), these lie at around 80% of the tip of the lobe. In the SR–SR lobe, the cut is taken at θ = −0.28. Here we see a similar behavior for the structure factor as we did in the θ = −0.15 cut of the lobe, but this time the superfluid fraction plays a much more important role. On the one hand, μ/J ≃ 2.4, both the superfluid fraction as well as the structure factor have sharp transitions that become even sharper at larger sizes. In fact, at T = 0 these transitions have been shown to be of first order, and the system goes directly from crystal to superfluid. Due to the finite temperature, they are continuous in our case. On the other hand, where μ/J ≃ 0.8, there appears a second-order phase transition into a supersolid region that spans all the way to μ/J = 0. Finally, we take a cut at θ = −0.23 of the LR–SR lobe. The behavior of this system seems to be quite different. The first thing to note is that the transitions on either side of the lobe are of second order. The other, and more important, observation is that now it appears that this system has supersolid behavior on both sides of the lobe: in addition to the expected supersolid at smaller μ/J, a region at μ/J above the crystal lobe appears where both the structure factor and the superfluid fraction are finite. If we recall figure 5, right panel, the LR–LR system showed that at θ = −0.15 as L increases the supersolid region disappears from the upper side of the lobe. In the case of the LR–SR system for θ = −0.23 the increase in the system size does not get rid of this supersolid phase.

**Figure 6.** SR–SR (left) at θ = −0.28 and LR–SR (right) at θ = −0.23. Solid blue: structure factor. Dashed black: superfluid fraction. Lines become thicker and darker as the system size goes up. In the SR–SR system, a supersolid at small μ/J persists at large systems, while at μ/J ≈ 2.4 the transition from crystal to superfluid becomes a direct first-order transition. For the LR–SR system, the curves change slowly at both sides of the crystal lobe, leading to persisting supersolids.
Download figure:
Standard image

3.3. Finite-temperature results

As a final calculation, we take a look at the important role that the temperature plays in both the melting of the crystal and the supersolid region. In this section, we will use the same cuts as in the previous section (θ = −0.28 for SR–SR, θ = −0.23 for LR–SR and θ = −0.15 for LR–LR) so that each system will posses all the possible phases: crystal, superfluid and supersolid. Each cut is investigated for 0.05 < T < 0.3 at a system size of L = 6. First, we analyze the structure factor to study how the ρ = 2/3 crystal melts with an increase in temperature (the second row of figure 7). For the SR–SR interactions, even at a temperature of 0.3 there still exists a bump in the structure factor which indicates that the crystal has not yet completely melted, while for both the long-ranged lobes the crystal melts by T ≃ 0.3. Interestingly, the SR–SR crystal and the LR–LR crystal are approximately the same size at T = 0.05, but by T = 0.3, one has melted while the other still exists. That means that the system with short-ranged interactions holds its crystal structure better at higher temperatures than does our system with all long-ranged interactions. Looking at the LR–SR lobe, we see that its crystal at this cut starts off smaller, yet it melts at about the same temperature as the one for LR–LR interactions. It seems that the dipolar repulsion helps stabilize the crystal structure over a larger temperature range, while the long-ranged hopping destroys the crystal more quickly because of the extra kinetic energy.

Maybe more importantly, we now study the melting of the supersolid for these same cuts. Figure 7 shows the structure factor, superfluidity and supersolidity as a function of temperature for each of the different systems. Since the supersolid is defined by having both a non-zero structure factor and non-zero superfluid fraction, by combining the graphs we are able to see where these regions exist and also how they melt with increased temperature (the bottom row of figure 7 shows a product of the structure factor and the superfluid fraction, which remains finite only where the two coexist). A common feature of all the graphs is the spikes on either side of the plateaus. These are regions where a phase transition occurs but it does not necessarily imply that a supersolid region exists. Most likely, these features appear due to the finite size of the system and the resulting smooth transitions of the superfluid fraction and structure factor. At larger sizes, the transitions would be much sharper at these points, the regions where a finite structure factor and superfluid fraction coexist would shrink and the spikes would diminish. From the previous section, we can assume that for the SR–SR system they would disappear completely at the upper transition from the crystal lobe while for the other two systems there would still exist a small supersolid region.

Returning to the main focus, the small-μ region, we see that in each case a supersolid region appears that extends from the left side of the plateau all the way to μ/J = 0. In every system, this supersolid region also exists for a finite range of temperatures. For SR–SR interactions, it gradually decreases but still extends all the way out past T = 0.3. For the LR–SR interactions, the supersolid region again slowly melts but now disappears at T ≃ 0.23, just below the spot where the crystal melted. The supersolid region for the LR–LR system appears to have the largest magnitude of the three systems, but then rapidly melts at T ≃ 0.3.

In order to compare these transitions more quantitatively, we take a cut along μ/J = 0 for each system, shown in figure 8. All three systems show a relatively similar and steady value for the structure factor. Hence, the values of the superfluid fraction are going to determine the existence of the supersolid regions. The first plot shows the SR–SR system at the θ = −0.28 cut, and we can see that the superfluid fraction stays non-zero all the way out to T = 0.35. The LR–SR system has a very similar behavior at the θ = −0.23 cut, but in this case the supersolid is nearly completely melted by T = 0.35. The final plot is the LR–LR system at θ = −0.15, which behaves slightly differently. The most important difference is that the starting value of the superfluid fraction is higher than in the first two plots. This should therefore make the supersolid region more pronounced. But even though the superfluid fraction has the highest value for this system, it decays more quickly and reaches values similar to the SR–SR system at T ≈ 0.35.

4. Conclusion

In this paper, we have presented a QMC study of dipolar spin models that describe various systems of ultracold atoms, molecules and ions. We have presented predictions concerning the phase diagram of the considered systems at zero and finite temperatures, and described the appearance and some properties of the superfluid, supersolid and crystalline phases. While the results are not surprising and resemble earlier obtained results for similar systems in 1D and 2D, the main advantage of our study is that it is directly relevant to the current experiments:

The results for the XXZ model with SR tunneling also apply to ultracold gases of polar bosonic molecules in the limit of hard-core bosons. Note that the earlier works [15, 16] have concentrated on the appearance of the supersolid phase and devil's staircase of crystalline phases in the square lattice [15], and supersolid and emulsion phases in the triangular lattice [16]. Here we focus on the hard-core-spin limit, and we compare it and stress differences with other models, such as those with LR tunneling, i.e. LR XX interactions.
The results for the XXZ models with LR tunneling apply to systems of trapped ions in triangular lattices of microtraps. These results are novel, since so far such models have been studied only using various techniques in 1D and using the mean-field approach in 2D. While the first experimental demonstrations of such models were restricted to a few ions (see, for instance, [3]), many experimental groups are working on an extension of such ionic quantum simulators to systems of many ions in microtraps [37]. In fact, very recently, the NIST group has engineered 2D Ising interactions in a trapped-ion quantum simulator with hundreds of spins [38]. Although in this experiment the quantum regime has not yet been achieved, it clearly opens the way toward quantum simulators of spin models with LR interactions. We expect that in the near future the result of our present theoretical study will become directly relevant for experiments.

Acknowledgments

This work was supported by the International PhD Projects Programme of the Foundation for Polish Science within the European Regional Development Fund of the European Union through agreement no. MPD/2009/6. We acknowledge financial support from the Spanish Government grants TOQATA (FIS2008-01236) and Consolider Ingenio QOIT (CDS2006-0019), EU IP AQUTE, EU STREP NAMEQUAM, ERC Advanced Grant QUAGATUA, CatalunyaCaixa, Alexander von Humboldt Foundation and Hamburg Theory Award. MM and JZ thank Lluis Torner, Susana Horváth and all ICFO personnel for hospitality. JZ acknowledges support from the Polish National Center for Science through project no. DEC-2012/04/A/ST2/00088.

Quantum spin models with long-range interactions and tunnelings: a quantum Monte Carlo study

Article metrics

Author e-mails

Author affiliations

Author notes

Dates

Abstract

1. Introduction

2. Review of the results in one dimension

3. Dipolar hard-core bosons on the triangular lattice

3.1. Vanishing tunneling

3.2. Low-temperature phase diagram at finite tunneling

3.3. Finite-temperature results

4. Conclusion

Acknowledgments

Footnotes

Quantum spin models with long-range interactions and tunnelings: a quantum Monte Carlo study

Article metrics

Share this article

Author e-mails

Author affiliations

Author notes

Dates

Abstract

1. Introduction

2. Review of the results in one dimension

3. Dipolar hard-core bosons on the triangular lattice

3.1. Vanishing tunneling

3.2. Low-temperature phase diagram at finite tunneling

3.3. Finite-temperature results

4. Conclusion

Acknowledgments

Footnotes