Evidence of spin liquid with hard-core bosons in a square lattice

We show that laser assisted hopping of hard core bosons in a square optical lattice can be described by an antiferromagnetic $J_{1}$-$J_{2}$ XY model with tunable ratio of $J_{2}/J_{1}$. We numerically investigate the phase diagram of the $J_{1}$-$J_{2}$ XY model using both the tensor network algorithm for infinite systems and the exact diagonalization for small clusters and find strong evidence that in the intermediate region around $% J_{2}/J_{1}\sim 0.5$, there is a spin liquid phase with vanishing magnetization and valence bond orders, which interconnects the Neel state on the $J_{2}\ll J_{1}$ side and the stripe antiferromagnetic phase on the $% J_{2}\gg J_{1}$ side. This finding opens up the possibility of studying the exotic spin liquid phase in a realistic experimental system using ultracold atoms in an optical lattice.

We show that laser assisted hopping of hard core bosons in a square optical lattice can be described by an antiferromagnetic J1-J2 XY model with tunable ratio of J2/J1. We numerically investigate the phase diagram of the J1-J2 XY model using both the tensor network algorithm for infinite systems and the exact diagonalization for small clusters and find strong evidence that in the intermediate region around J2/J1 ∼ 0.5, there is a spin liquid phase with vanishing magnetization and valence bond orders, which interconnects the Neel state on the J2 ≪ J1 side and the stripe antiferromagnetic phase on the J2 ≫ J1 side. This finding opens up the possibility of studying the exotic spin liquid phase in a realistic experimental system using ultracold atoms in an optical lattice.

PACS numbers:
A spin liquid phase is an exotic state of matter that does not break any symmetry of the Hamiltonian and has no conventional order even at zero temperature [1]. A number of microscopic Hamiltonians with frustrated quantum magnetic interaction could support a spin liquid phase [1][2][3][4][5][6]. In particular, very recently, numerical investigations based on complementary methods have found strong evidence that the antiferromagnetic J 1 -J 2 Heisenberg model may have a spin liquid phase in a square lattice [4,5]. On the experimental side, several materials are suspected to be in a spin liquid phase at very low temperature [1]. However, due to complication of physics in these materials, it is hard to make a direct connection of the prediction from the simplified microscopic models and the phenomenology observed in real materials [1]. Ultracold atoms in an optical lattice provides a clean platform to realize microscopic models to allow for controlled comparison between theory and experiments [7,8]. Proposals have been made to implement the frustrated magnetic models in an optical lattice [9,10] and various required configurations of the optical lattices have been realized experimentally [10]. However, the direct magnetic Heisenberg coupling, which comes from the higher-order super-exchange interaction, is very weak under typical experimental conditions [9,11]. It is still very challenging to reach the extremely low temperature required to observe the ground state of the magnetic Heisenberg model in an optical lattice.
In this paper, we show strong evidence that a spin liquid phase can emerge in an antiferromagnetic J 1 -J 2 XY model in a square lattice. The calculations are based on two complementary methods: the recently developed tensor network algorithm applied directly to infinite systems [12,13] and the exact diagonalization of small clusters which is combined with the finite size scaling to infer the phase diagram [14]. Both methods suggest that in a small region around J 2 /J 1 ≈ 0.5, magnetization and valence bond solid orders all vanish, indicting a spin liquid phase as the ground state. Different from a Heisenberg model, a XY model can be realized with hard-core bosons in an optical lattice. Through control of the laser assisted hopping in a square lattice [15], we propose a scheme to implement the effective antiferromagnetic couplings for both the neighboring and the next neighboring sites with a tunable ratio of J 2 /J 1 . In this implementation, both J 2 and J 1 are determined by the hopping rates of the hard-core bosons in an optical lattice, which is much larger than the conventional super-exchange interaction for ultracold atoms in the Heisenberg model [9,11]. The large J 1 -J 2 couplings open up the possibility to experimentally realize this model and observe its spin liquid phase based on the state-of-the-art technology.
The J 1 -J 2 XY model is represented by the Hamiltonian (1) where X, Y represent the Pauli operators σ x and σ y , i, j and i, j denote respectively the neighboring and the next neighboring sites in a square lattice as shown in Fig.  1(a). To realize this model with hard core bosons, we consider ultracold atoms in different hyperfine spins |a and |b loaded into alternating square lattices A and B as shown in Fig. 1. This configuration can be experimentally realized with the spin-dependent lattice potential [16]. Atoms in spins |a (or |b ) freely tunnel in the lattice A (or B) with the hopping rate t, however, a direct hopping between the A, B lattices is forbidden due to the spin-dependent potential shift. Instead, the inter-lattice hopping is introduced by the laser induced Raman transition as shown in Fig. 1(a). We use three Raman beams, with wave-vectors k 1 , k 2 , and k 3 and Rabi frequencies Ω 1 , Ω 2 , and Ω 3 , respectively. The directions of the laser beams are shown in Fig. 1(b) with ∆k 12 = k 1 −k 2 = k ∆ŷ and ∆k 13 = k 1 − k 3 = k ∆x . The lase induced interlattice hopping rates for the neighboring sites are then for the hopping along the x, y directions, respectively. Assume Ω 3 = −Ω 2 and the Wannier function w (x i , y i ) symmetric along the x, y directions, we have t x = −t y = t ′ (we can always choose t ′ > 0 by setting an appropriate relative phase between Ω 1 and Ω 3 ). If the on-site atomic repulsion U satisfies U ≫ t, t ′ , we have the hard-core constraint with at most one boson per site. The hard-core bosons in this square lattice are then described by the Hamiltonian (2) The hard core bosons a i , b j satisfy the same commutators as the Pauli operators σ − i , σ − j , so with the mapping a i −→ σ − i and b j −→ σ − j for the odd numbers of rows, and a i −→ −σ − i and b j −→ −σ − j for the even numbers of rows, the Hamiltonian (2) is mapped to the J 1 -J 2 XY model in Eq. (1) with J 1 = t ′ /2 > 0 and J 2 = t/2 > 0. Apparently, the ratio J 2 /J 1 is tunable by changing the magnitude of the Rabi frequencies Ω * 1 Ω 3 . In the following, we calculate the phase diagram of the Hamiltonian (1) as a function of the dimensionless parameter J 2 /J 1 (J 1 is taken as the energy unit). In the limit J 2 /J 1 ≪ 1, the J 1 term dominates and the ground state is magnetized with a Neel order at the momentum k = (π, π). In the opposite limit J 2 /J 1 ≫ 1, the ground state has a stripe magnetic order at the momentum (π, 0) or (0, π), which minimizes the energy of the J 2 term. In the intermediate region with J 2 /J 1 ∼ 0.5, the Hamiltonian is highly frustrated with competing interaction terms. Our main purpose is to find out the phase diagram in this region through controlled numerical simulations.
Our numerical simulations are based on two complimentary methods: exact diagonalization (ED) for small clusters [14] and tensor network simulation for infinite systems [12,13]. The ED method is limited by the cluster size, and we use extrapolation based on the finite-size scaling to infer the phase diagram for the infinite system. The tensor network algorithm is an recently developed simulation method inspired by quantum information theory [12]. It can be considered as an extension of the density matrix renormalization group (DMRG) method to the two dimensional case, replacing the matrix product state in the DMRG method with the tensor network state that better matches the geometry of the underlying lattice [12]. We use a particular version of the tensor network algorithms, the infinite projected entangled pair states (iPEPS) method [13], which applies directly to infinite systems using the translational symmetry. To take into account the ordered states for the Hamiltonian (1) that spontaneously break the translational symmetry, in our simulation we take a unit cell (typically 2 × 2 and 4 × 4) that is large enough to incorporate the relevant symmetry breaking orders [17]. We apply imaginary time evolution to reach the ground state of the Hamiltonian. To avoid being stuck in a metastable state, we take a number of random initial states for the imaginary time evolution and pick up the ground state as the one which has the minimum energy over all the trials. The accuracy of the iPEPS simulation depends on the internal dimension D of the tensor network state. The simulation time scales up very rapidly with the dimension D, which limits D to a small value in practice. We typically take D between 4 to 6 in our simulation. Figure 2 shows the major result from the iPEPS simulation. First, we look at the average magnetization where the average is taken over the N s sites in the unit cell. The calculation shows that for small or large J 2 /J 1 , the ground states are magnetic (with the Neel or the stripe order, respectively), which is consistent with our intuitive picture. In the intermediate region with 0.46 J 2 /J 1 0.54, there is a sudden drop of all the magnetic orders to a tiny value. Although the iPEPS method under a small dimension D could be biased toward a less entangled state, which is typically an ordered state, it would not be baised toward a disordered spin liquid state. So, when we see a big sudden drop of the magnetic orders from the simulation, it must be a real effect, strongly indicating there is a new phase in the intermediate region with vanishing magnetic orders. The remaining small m s may be due to the finite dimension D and should vanish when D is scaled up.
To figure out the property of the phase in the intermediate region, we further check different kinds of valence bond solid orders. We calculate all the neighboring valence bonds σ i · σ j in the unit cell and the result is shown in Fig. 2. For a valence bond solid state, the spatial symmetry should be spontaneously broken for the valence-bond distribution. Figure 2 shows that in the entire region of J 2 /J 1 , the valence bond distribution has the same symmetry as the underlying Hamiltonian, which indicates that the ground state of the Hamiltonian (1) has no valence bond solid orders. Together with the above calculation of the magnetic orders, this suggests that the  [4,5].
To further confirm this picture, we calculate the longrange spin correlation and dimer correlation with the iPEPS method and the result is shown in Fig. 3 for J 2 /J 1 = 0.1, 0.5 and 0.9. The spin correlation σ i · σ j is calculated along the diagonal directionx +ŷ . Both the Néel and the stripe phases have long-range correlations, with constant or staggered values along the diagonal direction. The intermediate phase has an exponentially decaying spin-spin correlation, which is in agreement with the behavior of the Z2 spin liquid phase with a finite spin gap [1,4]. The dimer operator D α i is defined by D α i = σ i · σ i+α for the bond (i, i + α), where α =x orŷ denote the orientation of the dimer. In Fig.  3(c), we show the dimer-dimer correlations ∆D x i ∆D x j and ∆D y i ∆D y j at J 2 /J 1 = 0.5 along the diagonal direction. The correlations are exponentially decaying with distance, in agreement with a spin liquid phase with no dimer orders.
In the following, we present study of the Hamiltonian (1) with the complementary ED method, which provides further evidence for a spin liquid phase in the intermediate region. To be consistent with the periodic boundary condition required for the finite size scaling and to incorporate the momentum k = (π, π) responsible for the Neel order, the size of the clusters for the ED is taken to 16, 20 and 32 sites. From the spin correlation σ i · σ j , we calculate the corresponding static structure factor m 2 s (k, N ) = (1/N ) ij e ik·(ri−rj ) ∆σ i · ∆σ j , where N is the size of the cluster and ∆σ i ≡ σ i − σ i . The Neel order and the stripe order correspond to peaks at k = (π, π) and (π, 0), respectively. Finite-size clusters always have non-zero order parameters, and one needs to do finite size scaling, with a simple scaling formula m 2 s (k, N ) = m 2 s (k, ∞) + a/ √ N ( √ N corresponds to the linear size), to infer the value of m 2 s (k, ∞) for the infinite system. In Fig. 4, we show the finite size scaling for m 2 s (k, N ) at J 2 /J 1 = 0, 0.5, and 0.9 in three different regions. The results are consistent with the findings from iPEPS method, i.e., there is a stripe order with k = (π, 0) at J 2 /J 1 = 0.9 and a Neel order with k = (π, π) at J 2 /J 1 = 0. At J 2 /J 1 = 0.5, the finitesize scaling indicates a vanishing stripe order. However, at k = (π, π), the data become non-monotonic with N due to the shape of the cluster and the finite-size scaling becomes inconclusive in this case. The non-monotonic shape effect has also been observed in ED for the J 1 -J 2 Heisenberg model [14].
To check for possible valence bond solid orders from ED, we similarly calculate the structure fac- clockwise (anticlockwise) cyclic permutation operator on the plaquette i with its explicit (lengthy) expression given in [17,18]. The rotational symmetry is always preserved at finite size, so we only need to check one component of the dimer order, say D x i . At finite size, the structure factors peak at k = (π, 0) for the dimer order D x i and at k = (π, π) for the plaquette order P i , however, an extrapolation to the infinite system at these momenta as shown in Fig. 5 indicates vanishing dimer and plaquette orders in all three regions of J 2 /J 1 . This result, again, is in agreement with the finding from the iPEPS calculation.
Before concluding the paper, we briefly discuss the experimental signature of the three different phases for the Hamiltonian (1) in the implementation with hard-core bosons. The Neel ordered state and the strip phase cor-respond to Bose-Einstein condensates at the momenta k = (π, π) and k = (π, 0), respectively. The standard time-of-flight imaging measurement can then reveal the condensate peak at these nontrivial momentum points [8]. The spin liquid phase, on the other hand, would not show any condensation peaks due to lack of magnetic orders. Furthermore, it has a spin gap which implies a charge gap in implementation with hard-core bosons. We therefore expect to see an incompressible phase at half filling, which is different from the Mott insulator state at the integer filling. It is also be distinguished from a charge density wave state since the density distribution in this case is still homogeneous without any solid order.
In summary, we have proposed an experimentally feasible scheme to implement the J 1 -J 2 XY model with ultracold hard-core bosons in a square optical lattice. Through detailed numerical simulation of this model using two complementary methods, we find strong evidence that this model has a spin liquid phase in the intermediate region of J 2 /J 1 . The proposed experimental implementation, with tunable ratio of J 2 /J 1 , opens up a realistic possibility to look for the long-pursued spin liquid phase in a well controlled Hamiltonian model.