# Exotic quantum spin models in spin-orbit-coupled Mott insulators

## Abstract

We study cold atoms in an optical lattice with synthetic spin-orbit coupling in the Mott-insulator regime. We calculate the parameters of the corresponding tight-binding model using Peierls substitution and “localized Wannier states method” and derive the low-energy spin Hamiltonian for fermions and bosons. The spin Hamiltonian is a combination of Heisenberg model, quantum compass model and Dzyaloshinskii-Moriya interactions and it has a rich classical phase diagram with collinear, spiral and vortex phases.

Since the first experimental realization of Bose-Einstein condensate (BEC), cold atoms have proven to be an excellent playground for studying many-body physics Lewenstein et al. (); Bloch et al. () and many interesting phenomena take place for strongly interacting atoms in an optical lattice. These studies started with an experimental observation of superfluid to Mott-insulator phase transition Greiner et al. () which was followed by experimental and theoretical work on both Bose and Fermi gases in lattices of different dimensionality and in various parameter regimes Lewenstein et al. (); Bloch et al. (). The key features of cold atoms in optical lattices are the excellent tunability of parameters and the fact that the sample is almost perfectly described by the Hubbard model in the deep lattice regime Bloch et al. (). Since it is well known that the Hubbard model is mapped to an effective spin Hamiltonian for Mott-insulator phases with integer filling MacDonald et al. (), it is clear that cold atoms can be used to “engineer” various quantum spin systems Kuklov and Svistunov (); Duan et al. () from those described by the Heisenberg model to more exotic ones, like the Kitaev model Kitaev (). In designing effective spin systems different tools like polar molecules Micheli et al. () and tilted optical lattices Simon et al. () can also be used. In recent years there has been a lot of interest in creating artificial Abelian and non-Abelian gauge fields in cold-atom systems Dalibard et al. () and successful experimental realizations of synthetic magnetic Lin et al. (a); Aidelsburger et al. () and electric field Lin et al. (b) and spin-orbit coupling (SOC) Lin et al. (c) have been reported. The role of SOC in cold atoms has been extensively studied following the theoretical proposal for the creation of artificial SOC Stanescu et al. (a); Ruseckas et al. () and rich phase diagrams have been found in BECs Stanescu et al. (b); Wang et al. (); Sinha et al. () and fermionic systems Sato et al. (); Iskin and Subasi (); Takei et al. ().

In this letter we combine optical lattice and SOC which, in the deep lattice regime, leads to tight-binding description with non-zero “spin-flip” hopping between neighboring sites Goldman et al. (a). We show that in the Mott-insulator phase with integer filling the system is described by an interesting effective spin Hamiltonian which is a combination of Heisenberg model, quantum compass model and Dzyaloshinskii-Moriya terms. We note that combination of an optical lattice and SOC has already been considered with a purpose of studying superfluid-insulator transition Graß et al. (), topological phase transitions Goldman et al. (b) and BEC dynamics Larson et al. (). In the context of solid-state physics spin models resulting from Mott-insulators with strong SOC were studied in Ref. Stekhtman et al. (a, b); Zheludev et al. (); Jackeli and Khaliullin ().

We study a two-dimensional system of pseudospin- atoms on a square optical lattice with synthetic SOC. The single-particle physics is described by the Hamiltonian:

(1) |

with and being the atomic momentum and mass; , the lattice depth in and direction, ( is the lattice spacing), the unit matrix, Pauli matrices; and characterize the SOC. Since (1) is invariant under lattice translations, its eigenstates have Bloch-wave form: and , being a lattice vector. We are interested in the deep lattice regime in which pairs of energy bands are well separated and the low-energy physics is captured by the lowest pair of bands which touch at , , and in the energy spectrum. In this regime the system is well described by the tight-binding approximation in which the Hilbert space is spanned by states localized on individual lattice sites and the tunneling exists only between nearest-neighbor sites. There are two localized states per site (, ), hence we have two effective particle species. This is the most general tight-binding description of (1):

(2) |

where () creates (annihilates) a particle in the state , are the tunneling matrices and . In finding the elements of corresponding to , we use Peierls substitution and “localized Wannier states (LWS) method”. We write SOC in a gauge-field form: with and may use Peierls substitution to find tunneling matrices

(3) |

where are tunneling coefficients in the -direction in the absence of SOC, , ; and are dimensionless SOC strengths: , . However, Peierls substitution is only an approximation, valid for SOC weak with respect to the kinetic plus lattice part of . This is the case when , , where is the lattice recoil energy. While the SOC is quite weak in solid-state systems, in cold-atom systems it is typically very strong and in that case Peierls substitution is not completely valid. For example, if we combine optical lattice with spacing and Rashba SOC generated by a scheme described in Ref. Campbell et al. (), we obtain , while the usual experimental values of lattice depth are (for smaller values of the tight-binding approximation is not valid). For the SOC scheme experimentally realized Lin et al. (c) . Since the validity condition for Peierls substitution is not completely satisfied, we calculate tunneling matrices using LWS method which is more involved and requires numerical approach but it does not contain any approximation. In a system with sites and periodic boundary conditions, Wannier states for a single band are defined as Wannier (); Kohn (): , where are Bloch states and is an arbitrary phase. In the absence of SOC it is possible to obtain maximally LWS by varying the phase of each Bloch state Kohn () and in the deep lattice regime these maximally LWS are well localized on individual sites. In the presence of SOC it is generally not possible to have Wannier states localized on individual sites if these are constructed from Bloch states of a single band. Therefore we consider generalized Wannier states introduced in Ref. Marzari and Vanderbilt (): , where are unitary matrices which mix Bloch states of the two bands. We find maximally LWS by minimizing the functional with respect to matrix elements of ( is an expectation value associated to ). The minimization is done numerically and example of an algorithm is given in Ref. Marzari and Vanderbilt (). Numerical results show that the tunneling matrices still have the form given in (3) but now the parameters , , and are some more general functions of , , , , and . It can be shown that the structure of tunneling matrices (3) follows from the symmetries of .

Peierls substitution has the advantage to give an analytical form for tunneling matrices, however it does not give any information about the Wannier states, whereas the LWS method explicitly gives states , which is important in interpreting the experimental data.

In Fig. 1 we compare tunneling amplitudes in the Rashba-coupling case obtained by Peierls substitution and LWS method for . They show excellent accord for small and sizable differences for larger ones.

Cold atoms in optical lattices are described by the tunneling Hamiltonian (2) plus interactions Bloch et al. ():

(4) |

where is a number of particles in state , are interaction coefficients and :(…): denotes normal ordering of creation and annihilation operators.

We are interested in the Mott-insulator regime with , and integer number of atoms per site ( for fermions and any integer for bosons). In this case interactions (4) are the dominant part in the full Hamiltonian and we may treat the problem perturbatively by taking as starting Hamiltonian and as perturbation. The ground state of is a state with uniform distribution of atoms and the ground state degeneracy is very large since there are two states per site that atoms can occupy. The perturbation couples the ground-state manifold and excited states of and the resulting low-energy effective Hamiltonian can be calculated in the second order of perturbation theory Kuklov and Svistunov (); Duan et al. ():

(5) |

where and label states in the ground-state manifold, while labels the excited states of .

We first calculate the effective low-energy Hamiltonian for fermions. Since two fermions cannot occupy the same quantum state the only interesting regime is when . The only relevant interaction coefficient is and the excited states of are those with two fermions of different species at the same site. Now it is convenient to introduce isospin operators: , and . Using (5) we obtain

(6) |

where , , , , and are introduced in (3). The Hamiltonian is a combination of Heisenberg model, compass model and Dzyaloshinskii-Moriya-type terms; for (no SOC) and we recover the Heisenberg model MacDonald et al. (); Kuklov and Svistunov (); Duan et al. ().

Next we find the effective low-energy Hamiltonian for bosons, and now the number of atoms per site can be greater than one. The calculation for any and general interaction coefficients is very cumbersome, however it simplifies for , where or when Kuklov and Svistunov (). We express the Hamiltonian in terms of spin- operators defined in the same way as in the previous case. For , the Hamiltonian in the first order of is where and .

For and general we obtain

(7) |

where , , , , , and , are given in (6) with replaced by . Since the atoms usually used in experiments have almost spin-independent interactions we assume (): this simplifies the Hamiltonian and yields .

We intend to find the classical zero-temperature phase diagram of (7) with . Some previous papers presented models combining Heisenberg, Dzyaloshinskii-Moriya and compass-model interactions Stekhtman et al. (a, b); Zheludev et al. () but they did not provide a complete phase diagram, neither at a classical level, usually considering only small SOC. In our approach we treat the spins as classical vectors and we aim to find the configurations which minimize the energy, with constant . We did our computations usually on -site lattices and finite-size effects are negligible. In Fig. 2 we show the phases and in Fig. 3 the corresponding phase diagram. We obtain two Ising-type phases [ferromagnet (Fig. 2a) and stripes (Fig. 2e)], coplanar spirals (Fig. 2b) and three-dimensional ordered phases with vortices (Fig. 2c) or antivortices (Fig. 2d). In describing our results, it is helpful to focus on a so called “basic region” given by the triangle : the solutions for other parameters can be obtained by simple mappings, e.g. ground-state configurations in region are obtained by simultaneous -rotation of spins and sites of ground states in the “basic region”. Upon activating SOC, the ferromagnet is immediately replaced by spiral waves, whose spatial periodicity reduces from several to three sites upon increasing and ; we found both commensurate and incommensurate waves. When the compass-model term becomes dominant over the Dzyaloshinskii-Moriya one, another coplanar phase appears, the ferromagnetic stripe order, either directly or via the three-dimensional ordered phases. We always find non degenerate classical ground states, except along the dashed lines (Fig. 3) which indicate points in the parameter space with a continuous degeneracy of classical ground states. However, we expect this degeneracy to be removed by slight deviations of the realistic engineered SOC with respect to the Rashba-Dresselhaus form of the coupling Campbell et al. (). The dashed lines also represent the boundaries between phases with stripes of different orientation, i.e. between the phase shown in Fig. 2e and the one obtained by rotating the sites and spins of the latter by around the -axis. The vortex phase (Fig. 2c) takes place along the diagonal : vortices are left-handed in the region with smaller SOC and right-handed in the one with larger SOC. The antivortex phase (Fig. 2d) is found along the diagonal and the configuration (d) is obtained from the phase (c) by a transformation which reflects sites (but not spins) with respect to the -axis. For a better identification of the phase properties, we consider their behavior with respect to the breaking of the translational symmetry of (7). While all the phases (except the ferromagnet) break this symmetry, they do it in a different way, i.e. the stripe phase in Fig.2e is not invariant under one-lattice-site translation along -direction, but it is invariant under two-lattice-sites translations in -direction and under one-lattice-site translation in -direction; the phases with vortices or antivortices are not invariant under one- and two-lattice-sites translations in and -direction, but they are invariant under three-lattice-sites translations. Then we can understand the evolution from stripe to vortex phase as a transition in which two-lattice-sites translational symmetry becomes broken. The same reasoning applies for the rest of the phase diagram.

It is important to emphasize that the classical analysis yields no gapless modes in the whole parameter region except for the diagonal lines (, ) in the stripe region. The absence of these gapless modes provides stronger guidelines for further analysis in a semiclassical or quantum approach.

In summary, we studied the effects of spin-orbit coupling in cold atoms in an optical lattice in the Mott-insulating regime. We derived the tight-binding model using Peierls substitution and Localized Wannier State method and obtained the effective low-energy Hamiltonian for fermions and bosons: this takes the form of an exotic spin model with Heisenberg, compass-model and Dzyaloshinskii-Moriya interactions. We determined the classical phase diagram for this model and showed that the interplay between the different interactions is responsible for a large variety of phases: ferromagnet, spirals, stripes, three-dimensional vortex and antivortex phases. We expect that our classification of ground states could generally survive in a quantum approach; in fact, except for some particular cases we mentioned in the discussion, on the classical level there are no degeneracies which would be lifted by quantum fluctuations.

We thank W.S. Cole, S. Zhang, A. Paramekanti and N. Trivedi for discussions. J.R. thanks Stephen Powell and Qinqin Lu for discussions. This research was supported by US-ARO, JQI (A.D.C. and V.G.), ARO-MURI (J.R.) and JQI-NSF-PFC (K.S.).

### References

- M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. S. De, and U. Sen, Adv. Phys. 56, 243 (2007).
- I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
- M. Greiner, M. O. Mandel, T. Esslinger, T. Hänsch, and I. Bloch, Nature 415, 39 (2002).
- A. H. MacDonald, S. M. Girvin, and D. Yoshioka, Phys. Rev. B 37, 9753 (1988).
- A. B. Kuklov and B. V. Svistunov, Phys. Rev. Lett. 90, 100401 (2003).
- L.-M. Duan, E. Demler, and M. D. Lukin, Phys. Rev. Lett. 91, 090402 (2003).
- A. Y. Kitaev, Ann. Phys. 321, 2 (2006).
- A. Micheli, G. K. Brennen, and P. Zoller, Nat. Phys. 2, 341 (2006).
- J. Simon, W. S. Bakr, R. Ma, M. E. Tai, P. M. Preiss, and M. Greiner, Nature 472, 307 (2011).
- J. Dalibard, F. Gerbier, G. Juzelinas, and P. Öhberg, Rev. Mod. Phys. 83, 1523 (2011).
- Y.-J. Lin, R. L. Compton, K. Jimenez-Garcia, J. V. Porto, and I. B. Spielman, Nature 462, 628-632 (2009).
- M. Aidelsburger, M. Atala, S. Nascimbene, S. Trotzky, Y.-A. Chen, and I. Bloch, Phys. Rev. Lett.107, 255301 (2011).
- Y.-J. Lin, R. L. Compton, K. Jimenez-Garcia, W. D. Phillips, J. V. Porto, and I. B. Spielman, Nature Physics 7, 531-534 (2011).
- Y.-J. Lin, K. Jimenez-Garcia, and I. B. Spielman, Nature 471, 83-86 (2011).
- T. D. Stanescu, C. Zhang, and V. M. Galitski, Phys. Rev. Lett. 99, 110403 (2007).
- J. Ruseckas, G. Juzelinas, P. Öhberg, and M. Fleischhauer, Phys. Rev. Lett. 95, 010404 (2005).
- T. D. Stanescu, B. Anderson, and V. Galitski, Phys. Rev. A 78, 023616 (2008).
- C. Wang, C. Gao, C.-M. Jian, and H. Zhai, Phys. Rev. Lett. 105, 160403 (2010).
- S. Sinha, R. Nath, and L. Santos, Phys. Rev. Lett. 107, 270401 (2011).
- M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. Lett. 103, 020401 (2009).
- M. Iskin and A. L. Subasi, Phys. Rev. A 84, 043621 (2011).
- S. Takei, C.-H. Lin, B. M. Anderson, and V. Galitski, Phys. Rev. A 85, 023626 (2012).
- N. Goldman, A. Kubasiak, A. Bermudez, P. Gaspard, M. Lewenstein, and M. A. Martin-Delgado, Phys. Rev. Lett. 103, 035301 (2009).
- T. Graß, K. Saha, K. Sengupta, and M. Lewenstein, Phys. Rev. A 84, 053632 (2011).
- N. Goldman, W. Buegeling, and C. M. Smith, Europhys. Lett. 97, 23003 (2012).
- J. Larson, J.-P. Martikainen, A. Collin, and E. Sjöqvist, Phys. Rev. A 82, 043620 (2010).
- L. Stekhtman, O. Entin-Wohlman, and A. Aharony, Phys. Rev. Lett. 69, 836 (1992).
- L. Stekhtman, O. Entin-Wohlman, and A. Aharony, Phys. Rev. B 47, 174 (1993).
- A. Zheludev, S. Maslov, G. Shirane, I. Tsukada, T. Masuda, K. Uchinokura, I. Zaliznyak, R. Erwin, and L. P. Regnault, Phys. Rev. B 59, 11432 (1999).
- G. Jackeli and G. Khaliullin, Phys. Rev. Lett. 102, 017205 (2009).
- D. L. Campbell, G. Juzelinas, and I. B. Spielman, Phys. Rev. A 84, 025602 (2011).
- G. H. Wannier, Phys. Rev. 52, 191 (1937).
- W. Kohn, Phys. Rev. 115, 809 (1959).
- N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).