A continuum of compass spin models on the honeycomb lattice

Model Hamiltonians describing interacting spins localized on lattice sites are at the central stage in the field of quantum magnetism. A class of spin models, collectively known as the compass models [1], stand out owing to a unique feature they share in common: the spin exchange interactions differ for different lattice bond orientations. This is in contrast to the familiar Heisenberg model or the Ising model, where the exchange has the same form for all bonds connecting the nearest neighboring sites. The compass models arise naturally as low energy effective Hamiltonians in Mott insulators with orbital degrees of freedom [2–7] as well as interacting systems with spin–orbit coupling. These highly nontrivial models are also very appealing from a pure theoretical point of view because they offer a natural arena to study frustrated quantum magnetism [8, 9]. Exactly solvable compass models, the Kitaev model in particular, have played a pivotal role in stimulating the field of topological quantum computing [10, 11]. The rich physics contained in compass models has been reviewed recently in [1].

Our work is directly motivated by two well known compass models defined on the honeycomb lattice. The first example is the Kitaev model [11], where the exchange interactions between two neighboring sites are given by ${\sigma }_{i}^{x}{\sigma }_{j}^{x},{\sigma }_{i}^{y}{\sigma }_{j}^{y}$, and ${\sigma }_{i}^{z}{\sigma }_{j}^{z}$ respectively. As shown by Kitaev, this model is exactly solvable and has anyonic excitations obeying fractional statistics [11]. The spatially dependent exchange interactions suppress long-range spin order and support a quantum spin liquid (SL) ground state, one of the most sought after exotic many-body states in condensed matter physics [12]. The Kitaev model, despite its theoretical appeal, is neither readily realized in materials nor easily simulated with synthetical quantum matter such as cold atoms on optical lattices. Recently, the hybrid Kitaev–Heisenberg model, a linear superposition of a Kitaev term and a Heisenberg term, was proposed for iridium oxides and solved numerically [13]. Besides the SL phase, the hybrid model contains other interesting phases such as the stripe and the zigzag phase due to the competition between the two terms [13]. The phase diagram becomes even richer when off-diagonal spin exchange interactions are added [14–16].

The second example of compass models is the quantum 120^◦ model [6, 7]. It is very analogous to the Kitaev model but the spin operators ${\sigma }_{x},{\sigma }_{y}$, and ${\sigma }_{z}$ for the three bond directions are replaced by three (pseudo)spin 1/2 operators ${T}^{1},{T}^{2}$, and T³ respectively. They form an angle of 120° with each other on the xz plane in spin space, hence the name 'the 120^◦ model'. It was introduced to described the low energy physics of transition metal oxides [17] with doubly degenerate orbitals, e.g. orbital-only models of e_g orbitals on cubic lattice [3]. Later, two of us, and Wu independently, found that the 120^◦ model can be naturally realized in strongly interacting spinless p-orbital fermions on the honeycomb optical lattice [6, 7]. Although it is geometrically frustrated, spin wave analysis indicates that long-range order is stabilized by quantum fluctuations through the order by disorder mechanism [6, 7]. While the semiclassical spin-wave analysis is suggestive, the ground state of the 120^◦ model on honeycomb lattice remains to be identified unambiguously.

Given the apparent similarities between the Kitaev model and the 120^◦ model, it is natural to seek the conceptual and quantitative link between them. Indeed, these two models can be viewed as two instances of a more general class of compass models [18]. In this paper, we provide a concrete construction and propose a 'super model' which contains three paradigm models, the Ising, the Kitaev and the 120^◦ model, as special limits. It only has a single tuning parameter ${\theta }^{\prime }$ and a simple, intuitive picture for the three (pseudo)spin operators: they form a tripod in spin space as shown in figure 1. Analogous to tuning the tripod to raise or low a mounted camera in photography, dialing the angle ${\theta }^{\prime }$ between the three legs takes the Ising model (tripod fully closed) smoothly to the Kitaev model (tripod open with three legs orthogonal to each other) and then to the 120^◦ model (tripod fully open with three legs in the same plane). Immediately, one conjectures that the phase diagram of this continuum compass model is highly nontrivial containing drastically different long-range ordered states as well as SLs.

Zoom In Zoom Out Reset image size

We obtain the phase diagram of this 'super model' using tensor network states which have gained considerable success recently in the study of frustrated magnetism [19–22]. The results are summarized in figures 1 and 2. The order parameters are calculated using the tensor renormalization group (TRG) method formulated in thermodynamic limit [23, 24]. We show that the ground state of the quantum 120^◦ model is a long-range ordered dimer phase, and a SL phase exists in an extended region in our phase diagram⁵ . The numerical results of TRG are further confirmed and crosschecked with projected entangled pair states (PEPS) calculations [25, 26] for finite systems, exact diagonalization (ED), and spin wave analysis. We discuss the qualitative features of the quantum phase transitions between the SL phase and the dimer phase by introducing a topological charge (spin vortex) for the dimer configuration. We further show that the proposed tripod model can in principle be simulated with Hubbard model in the Mott insulating regime, e.g., using cold atoms on optical lattice with artificial gauge fields.

Zoom In Zoom Out Reset image size

We generalize the Kitaev and the ${120}^{\circ }$ model to the following continuum compass model defined on the two-dimensional honeycomb lattice

$$\begin{eqnarray}&&H(\theta )=J\displaystyle \sum _{{\bf{r}},\gamma }{T}_{{\bf{r}}}^{\gamma }(\theta ){T}_{{\bf{r}}+{{\bf{e}}}_{\gamma }}^{\gamma }(\theta )\end{eqnarray} \tag{ 1 }$$

where $J\gt 0$, and for each lattice site ${\bf{r}}$ the spin 1/2 operators are defined as

$$\begin{eqnarray}&&{T}^{\gamma }(\theta )=\displaystyle \frac{1}{2}({\tau }^{z}\mathrm{cos}{\phi }_{\gamma }+{\tau }^{x}\mathrm{sin}{\phi }_{\gamma })\mathrm{cos}\theta +\displaystyle \frac{1}{2}{\tau }^{y}\mathrm{sin}\theta \end{eqnarray} \tag{ 2 }$$

with the Pauli matrices ${\tau }^{x},{\tau }^{y},{\tau }^{z}$. Each site ${\bf{r}}$ is coupled to its neighbors ${\bf{r}}+{{\bf{e}}}_{\gamma }$, where ${{\bf{e}}}_{\gamma },\gamma =1,2,3$, denotes the three bond vectors of the honeycomb lattice. Geometrically, the three ${T}^{\gamma }$ form a tripod in the spin space as shown in figure 1: they are tilted from the xz plane by angle θ and, when projected onto the xz plane, are orientated along the corresponding bond direction ${{\bf{e}}}_{\gamma }$, i.e. at azimuthal angle ${\phi }_{\gamma }=0,2\pi /3,4\pi /3$ respectively. While ${T}^{\gamma }$ is most naturally defined through the tilting angle θ, it is much more convenient to introduce another angle, ${\theta }^{\prime }$, to discuss the various limits of $H(\theta )$. ${\theta }^{\prime }$ is the angle between T¹ and T², i.e., the two adjacent legs of the tripod. And it is related to θ by trigonometry

$$\begin{eqnarray}&&\mathrm{cos}{\theta }^{\prime }=1-\displaystyle \frac{3}{2}{\mathrm{cos}}^{2}\theta .\end{eqnarray} \tag{ 3 }$$

We will take ${\theta }^{\prime }$ as the only tuning parameter in the tripod model.

Three special limits of this model can now be identified. First of all, when ${\theta }^{\prime }=0,{T}^{\gamma }$ all collapse to ${\tau }_{y}$. The tripod is closed, and $H(\theta )$ is nothing but the Ising model

$$\begin{eqnarray}&&{H}_{{\rm{I}}}=\displaystyle \frac{J}{4}\displaystyle \sum _{{\bf{r}},\gamma }{\tau }_{{\bf{r}}}^{y}{\tau }_{{\bf{r}}+{{\bf{e}}}_{\gamma }}^{y}.\end{eqnarray} \tag{ 4 }$$

Note that we choose ${\tau }_{y}$ as the vertical axis in spin space instead of the usual convention of ${\tau }_{z}$ so that $H(\theta )$ reduces exactly to the ${120}^{\circ }$ model defined in our earlier work [6]. Secondly, when ${\theta }^{\prime }=90^\circ ,H(\theta )$ reduces to the Kitaev model since the operators ${T}^{\gamma }$ are now perpendicular to each other in the spin space. We can simply identify them as ${\sigma }^{x},{\sigma }^{y}$ and ${\sigma }^{z}$ (apart from a factor 1/2) in a rotated coordinate system as illustrated in figure 1(b). Thirdly, for ${\theta }^{\prime }=120^\circ ,H(\theta )$ becomes the quantum 120^◦ with ${T}^{1},{T}^{2},{T}^{3}$ all confined within the xz plane. It can be visualized as a fully open tripod.

As well known, the Ising model has antiferromangetic (AF) order with the order parameter

$$\begin{eqnarray}&&{O}_{1}=\langle | {\tau }^{y}| \rangle /2.\end{eqnarray} \tag{ 5 }$$

On the other hand side, the quantum 120^◦ model is conjectured to be long-range ordered despite the geometric frustration. We introduce the following 'order parameter' to measure the in-plane magnetization

$$\begin{eqnarray}&&{O}_{2}=\sqrt{\langle {\tau }^{x}{\rangle }^{2}+\langle {\tau }^{z}{\rangle }^{2}}/2.\end{eqnarray} \tag{ 6 }$$

By solving $H(\theta )$ using tensor network algorithms, we compute the average spin $\langle {\tau }_{{\bf{r}}}^{x,y,z}\rangle$ in the ground state. The main results are summarized in the schematic phase diagram in figure 1(b). Figure 2 shows the variation of the two order parameters introduced above as ${\theta }^{\prime }$ is changed. The region at small ${\theta }^{\prime }$ corresponds to the familiar Néel order which is characteristic of the classical Ising model and illustrated in the left inset of figure 2. Despite the increased quantum fluctuations as ${\theta }^{\prime }$ is increased, the Néel ordered phase persists up to ${\theta }^{\prime }\sim 87^\circ$. At the opposite end of large ${\theta }^{\prime }$, we find that the long-range spin order consists of a set of 'dimers,' i.e. opposite spins on neighboring sites, arranged into a periodic pattern of triangular lattice (figure 5(a)). The triangular lattice and its enlarged unit cell becomes transparent if we introduce a topological charge (red dot in figure 5(a)) for each hexagon with spins all pointing outwards. If we focus on one individual hexagon, e.g. the one shown in the right inset of figure 2, the orientations of the dimers happen to be also 60° (or equivalently 120°) apart. We will refer to this phase simply as the 'dimer phase.' In particular, it is the ground state of the quantum ${120}^{\circ }$ model on honeycomb lattice. This point will be further discussed in section 5.

Sandwiched between the Néel ordered phase and the dimer phase, a quantum SL phase is stabilized for ${\theta }^{\prime }\in [87^\circ ,94^\circ ]$. The conclusion is mainly based on the observation from figure 2 that the order parameters ${O}_{\mathrm{1,2}}$ in this region are nearly zero compared to those in other two phases. This conclusion is also consistent with the exactly solution of Kitaev model for ${\theta }^{\prime }=90^\circ$. The order parameters as functions of ${\theta }^{\prime }$ in figure 2 also suggest that the two quantum phase transitions in the tripod model may be qualitatively different. The gradual drop of O₂ at ${\theta }^{\prime }\gt 90^\circ$ indicates a continuous phase transition between the dimer phase and the SL phase. In contrast, the drop of the order parameter O₁ at ${\theta }^{\prime }\lt 90^\circ$ is rather sharp, pointing to a likely first-order phase transition. The details of the calculations leading to these results will be discussed below in section 3.

Recent developments of entanglement-based tensor network algorithms provide a novel, accurate approach to strongly correlated electron systems [25–29]. Particularly, they have been successfully applied to frustrated quantum magnets [19–22] and the t − J model [30, 31] to yield insights previously unattainable from conventional methods. To find the phase diagram of the proposed tripod model, we employ two complementary tensor networks algorithms, one for finite-size systems and the other for infinite systems in the thermodynamic limit, to find the ground state and the order parameters. In both algorithms, the ground state wave function is constructed as a network of local tensors defined on lattice sites. Each tensor has one physical index representing the spin degree of freedom and three virtual indices, each with bond dimension D, describing the quantum entanglement with its three neighboring sites.

We first apply the finite PEPS algorithm [25–27] to solve $H(\theta )$ for a six-site system with periodic boundary conditions. The ground state energies obtained coincide with those from ED. This suggests that PEPS is intrinsically superior compared to mean field theories when applied to frustrated spin Hamiltonians such as $H(\theta )$. The order parameters decay to zero as the ground state is approached for a finite system. Nonetheless, their decay behaviors are quite disparate for ${\theta }^{\prime }$ values in the Ising, Kitaev, and ${120}^{\circ }$ regions, suggesting three different phases. The details of the calculation are presented in appendix A.

To study the tripod model in the thermodynamic limit, we first find the converged ground state using imaginary time evolution and following the simple update scheme as described in [32] which generalizes the time-evolving block decimation (TEBD) [33] technique to two-dimensions. For a n-site unit cell, e.g. a six-site unit cell shown in figure 3(a), we need $3n/2$ different bond vectors that represent, roughly speaking, a mean-field approximation of the environment. Using these bond vectors, the simple update starts with n random tensors and iterates until convergence is achieved. At the end of the calculation, the ground state $| {\rm{\Psi }}\rangle$ is characterize by n tensors ${T}_{j},j\in [1,n]$. We then evaluate the expectation value of operator $O,\langle O\rangle =\langle {\rm{\Psi }}| O| {\rm{\Psi }}\rangle /\langle {\rm{\Psi }}| {\rm{\Psi }}\rangle$ which involves the (infinite) product of tensors T_j, using a real space coarse graining procedure known as higher-order TRG (HOTRG) [24] schematically shown in figure 3. We outline the main steps here. At the ith step of HOTRG, a local tensor, say Tⁱ₁, is regrouped with its three nearest neighbor tensors (${T}_{2},{T}_{4},{T}_{6}$) to form a new tensor ${\tilde{T}}_{1}^{i+1}$. More generally, for odd or even sites

$$\begin{eqnarray}&&{\tilde{T}}_{o}^{i+1}=\displaystyle \sum _{{\rm{s}}.{\rm{l}}.}{T}_{o}^{i}{T}_{2}^{i}{T}_{4}^{i}{T}_{6}^{i},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\tilde{T}}_{{\rm{e}}}^{{\rm{i}}+1}=\displaystyle \sum _{{\rm{s}}.{\rm{l}}.}{T}_{{\rm{e}}}^{{\rm{i}}}{T}_{1}^{i}{T}_{3}^{i}{T}_{5}^{i},\end{eqnarray} \tag{ 8 }$$

where the summation is over the shared legs (abbreviated as s.l., the solid lines in figure 3(b)) of the neighboring tensors, i.e. tensor contractions. The new tensors, each of which contains four old tensors, are of higher dimensions and truncated to have the same dimension as T_i via

$$\begin{eqnarray}&&{T}_{o/e}^{i+1}=\displaystyle \sum _{{\rm{s}}.{\rm{l}}.}{\tilde{T}}_{o/e}^{i+1}{U}_{x}{U}_{y}{U}_{z},\end{eqnarray} \tag{ 9 }$$

where the three projection tensors ${U}_{x,y,z}$, shown in figure 3(b), are obtained as follows. Take U_x as an example. First, ${\tilde{T}}_{j}$ is reshaped into matrix ${{ \mathcal T }}_{j}$ with the row corresponding to the leg along the x-direction. Then, a matrix ${ \mathcal U }$ is obtained by singular value decomposition

$$\begin{eqnarray}&&\displaystyle \sum _{j=1}^{6}{{ \mathcal T }}_{j}{{ \mathcal T }}_{j}^{\dagger }={ \mathcal U }{\rm{\Lambda }}{{ \mathcal U }}^{\dagger }.\end{eqnarray} \tag{ 10 }$$

Finally, U_x is obtained by truncating ${ \mathcal U }$ to a given truncation dimension χ and reshaping it back to the tensor form. The new tensors ${T}^{i+1}$ now form exactly the same honeycomb lattice structure as the old tensors Tⁱ but represent a larger system, see figure 3(c). This constitutes a single RG step. By iterating the RG steps many times, the converged result of $\langle O\rangle$ well approximates the expectation value in the thermodynamic limit.

Zoom In Zoom Out Reset image size

By following these procedures, we have calculated the ground state energy and the ground state expectation values of the order parameters ${O}_{1},{O}_{2}$ for different unit cell sizes, $n=2,4,6,8$. We found that the six-site unit cell gives the lowest energy. The two-site unit cell yields results in agreement with the six-site unit cell within the parameter region ${\theta }^{\prime }\lt 90^\circ$. The four-site and eight-site unit cells, however, lead to excited states with significantly higher energy. Thus, we conclude that the six-site unit cell is the most reasonable choices for all the ${\theta }^{\prime }$ values in the ground state calculation. In practice, one can safely use the two-site unit cell for ${\theta }^{\prime }\lt 90^\circ$ since it is significantly cheaper. The phase diagram and the spin configurations in the ordered phases shown in figure 2 are obtained by using the two-site unit cell for ${\theta }^{\prime }\leqslant 90^\circ$ and the six-site unit cell for ${\theta }^{\prime }\gt 90^\circ$.

To cross-check the TRG results, we perform the standard spin wave analysis of the tripod model. It is important to keep in mind that the validity of the spin wave theory, which can be viewed as expansion in series of $1/S$, becomes questionable in the limit of $S=1/2$. Yet the analysis offers a rough picture of the role played by geometric frustration and how the Néel order and dimer order get destroyed by the increased quantum fluctuations. As we will show below, the estimations of the two quantum critical points from the spin wave theory turn out to be in broad agreement with the phase digram predicted by the tensor network algorithms.

The analysis starts by partitioning the honeycomb lattice into the A and B sublattice and introducing ${S}_{{\bf{r}}}=\pm {T}_{{\bf{r}}}$ for all sites on the A (B) sublattice. Then the tripod Hamiltonian acquires a suggestive form

$$\begin{eqnarray}&&H(\theta )=\displaystyle \frac{J}{2}\displaystyle \sum _{{\bf{r}}\in A,\gamma }{[{S}_{{\bf{r}}}^{\gamma }(\theta )-{S}_{{\bf{r}}+{{\bf{e}}}_{\gamma }}^{\gamma }(\theta )]}^{2}+\mathrm{const}.\end{eqnarray} \tag{ 11 }$$

Here we have promoted the spin 1/2 operator ${\boldsymbol{\tau }}/2$ to general spin operator ${\bf{S}}$ with spin quantum number S. It follows that classical ground states are massively degenerate (except for the Ising limit). Any spin configurations with ${S}_{{\bf{r}}}^{\gamma }(\theta )={S}_{{\bf{r}}+{{\bf{e}}}_{\gamma }}^{\gamma }(\theta )$, i.e., the projection of ${\bf{S}}$ along the bond direction being the same for any two neighboring sites, will minimize the classical energy. This is a well known feature of compass models, see the review [1]. The special case of the classical ${120}^{\circ }$ model on honeycomb lattice was previously discussed in [7, 34]. We will confine our spin wave analysis to the simple case of spatially homogeneous spin configurations ${{\bf{S}}}_{{\bf{r}}}={{\bf{S}}}_{0}$ as done in [34]. The direction of ${{\bf{S}}}_{0}$ is characterized by its polar angle φ measured from ${\tau }_{y}$ and its azimuthal angle α of ${{\bf{S}}}_{0}$ measured from ${\tau }_{z}$ in ${\tau }_{x}$–${\tau }_{z}$ plane. The corresponding classical ground state energy per unit cell is

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{0}}{{S}^{2}J}=-\displaystyle \frac{3}{2}(2{\mathrm{sin}}^{2}\theta {\mathrm{cos}}^{2}\varphi +{\mathrm{cos}}^{2}\theta {\mathrm{sin}}^{2}\varphi ).\end{eqnarray} \tag{ 12 }$$

It is interesting to note that, coincidently, at the Kitaev point, ${\theta }^{\prime }=90^\circ$ which corresponds to $\theta ={\mathrm{cos}}^{-1}\sqrt{2/3},{E}_{0}$ is completely flat and does not depend on φ. For ${\theta }^{\prime }\lt 90^\circ ,{E}_{0}$ is minimized when $\varphi =0$ or π, corresponding to the two degenerate states with spin up or down in the Néel ordered phase. In contrast, for ${\theta }^{\prime }\gt 90^\circ ,{E}_{0}$ is minimized when $\varphi =\pi /2$, i.e., ${{\bf{S}}}_{0}$ lies within the ${\tau }_{x}$–${\tau }_{z}$ plane. Therefore, the mean field theory above predicts that the tripod model has a phase transition exactly at the Kitaev point.

Applying the Holstein–Primakoff transformation [35] to $H(\theta )$ and expanding the resulting Hamiltonian of bosons to order $1/S$, we compute the quantum fluctuation correction to the ground state energy for the two long-ranged ordered states respectively and find

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{1}}{{SJ}}=\displaystyle \frac{1}{N}\displaystyle \sum _{{\bf{k}},\lambda }| {\omega }_{\lambda }({\bf{k}})| -{\rm{\Delta }},\end{eqnarray} \tag{ 13 }$$

where ${\rm{\Delta }}=-{E}_{0}/{S}^{2}J,N$ is the number of sites within the A sublattice, the ${\bf{k}}$ summation is over the first Brillouin zone of the A sublattice, and ${\omega }_{\lambda }({\bf{k}})$ describe the spin wave dispersion for branch $\lambda =1,2$ and they are given by the eigenvalues of the matrix

$$\begin{eqnarray}\displaystyle \frac{1}{4}\left[\begin{array}{cccc}2{\rm{\Delta }} & {\beta }_{3}^{*} & 0 & {\beta }_{1}^{*}\\ {\beta }_{3} & 2{\rm{\Delta }} & {\beta }_{2}^{*} & 0\\ 0 & -{\beta }_{2} & -2{\rm{\Delta }} & -{\beta }_{3}^{*}\\ -{\beta }_{1} & 0 & -{\beta }_{3} & -2{\rm{\Delta }}\end{array}\right].\end{eqnarray} \tag{ 14 }$$

The expressions for ${\beta }_{1},{\beta }_{2}$ and ${\beta }_{3}$ are lengthy and tabulated in appendix B. In what follows, we will discuss the energy $E(\varphi ,\alpha )={E}_{0}+{E}_{1}$ separately for the two distinct cases: ${\theta }^{\prime }\lt 90^\circ$ and ${\theta }^{\prime }\gt 90^\circ$.

The results for $E(\varphi )$ from the spin wave analysis are plotted in the upper panel of figure 4 for several values of ${\theta }^{\prime }$ corresponding to the Néel ordered phase. One notices that the fluctuations do not change qualitatively the mean field ground state. $E(\varphi )$ reaches minima still at $\varphi =0$ or π for small ${\theta }^{\prime }$. However, as ${\theta }^{\prime }$ is increased, the energy $E(\varphi )$ becomes flatter. Eventually, as ${\theta }^{\prime }={\theta }_{{\rm{c}}}^{\prime }\sim 75.0^\circ$ (the top curve of figure 4), the energies for $\varphi =\pi /4,3\pi /4$ with proper choice of α become degenerate with those for $\varphi =0,\pi$. This signals the destabilization of the Néel order by quantum fluctuations. This occurs around ${\theta }_{{\rm{c}}}^{\prime }$, before the Kitaev point is approached. Note that in figure 4, only the region $\varphi \in [0,\pi /4]\cup [3\pi /4,\pi ]$ is shown. Outside this region (and also for ${\theta }^{\prime }\gt {\theta }_{{\rm{c}}}^{\prime }$), the lowest order spin wave theory based on the Ising-like antiferromagnetic order becomes ill defined.

Zoom In Zoom Out Reset image size

For ${\theta }^{\prime }\gt 90^\circ$, the classical ground state is continuously degenerate with $\varphi =\pi /2$ but arbitrary $\alpha \in [0,2\pi ]$. Quantum fluctuations lift the degeneracy and select a long-range ordered ground state via the 'order by disorder mechanism'. Such mechanism for the special case of ${\theta }^{\prime }=120^\circ$, i.e. the quantum 120^◦ model, has been discussed before in [6, 7, 34]. As shown in the lower panel of figure 4, the same physical picture continues to hold for the tripod model for ${\theta }^{\prime }\leqslant 120^\circ$: the energy E is minimized at ${\alpha }_{n}=n\pi /3$ with integer n. However, $E(\alpha )$ becomes increasingly flat as ${\theta }^{\prime }$ is decreased. At the critical point ${\theta }^{\prime }={\theta }_{{\rm{d}}}^{\prime }\sim 94.6^\circ$, additional minima of E appear at $\alpha ={\alpha }_{n}+\pi /6$. This indicates that the long-range order gets destroyed and replaced by a new phase for ${\theta }^{\prime }\leqslant {\theta }_{{\rm{d}}}^{\prime }$.

We emphasize that the simple version of spin wave theory outlined above is only intended to estimate the lower and upper critical points for the SL phase, ${\theta }_{{\rm{c}}}^{\prime }$ and ${\theta }_{{\rm{d}}}^{\prime }$. It can be further improved to properly treat general classical spin configurations. We will not do it here because the large S expansion by itself cannot unambiguously determine the order for our model of $S=1/2$ in the region ${\theta }^{\prime }\gt 90^\circ$.

Previous theoretical studies of the quantum ${120}^{\circ }$ model on the honeycomb lattice gave conflicted results. The spin wave analysis of [6] assumed a homogenous ground state ${{\bf{S}}}_{{\bf{r}}}={{\bf{S}}}_{0}$ and found quantum fluctuations prefer ${\alpha }_{n}=n\pi /3$. This led the authors to suggest that the ground state may be a simple ferromagnet of ${\bf{S}}$ with any choice of ${\alpha }_{n}$ (i.e. antiferromagnetic order in terms of the original spin ${\bf{T}}$ or ${\boldsymbol{\tau }}$). Reference [7] considered more general (inhomogeneous) classical ground states and discovered that, within spin wave theory, the ferromagnetic state is energetically less competitive than a 'fully packed unoriented loop configuration' with the same ${\alpha }_{n}$ values. The reason is quite subtle but argued to be physically robust: the loop configuration hosts maximum number of zero modes. This result obtained from semiclassical large-S expansion was conjectured to survive in the limit of $S=1/2$, i.e. the quantum ${120}^{\circ }$ model has a ground state with the six-site plaquette order [7]. However, no evidence of long-range order was found in ED studies where the spin correlation functions were computed for finite size clusters with periodic boundary conditions [34]. Instead, the ED results supported a trial wave function similar in spirit to the short-range resonating valence bond state, i.e., a liquid state with linear superposition of dimer covering of the lattice. Therefore, the true ground state of the quantum ${120}^{\circ }$ model was not settled.

Compared to these previous works, the numerical tensor network algorithm used here takes into account quantum fluctuations beyond the lowest order spin wave theory, works directly in the thermodynamic limit, and starts with unbiased (random) choice of tensors as variational parameters. It is capable of describing both the long-range ordered and the spin-disordered ground states. We find the ground state of the quantum ${120}^{\circ }$ model is the dimer phase illustrated in figure 5(a) where the arrows denote the direction $\hat{n}$ of spin average $\langle {\boldsymbol{\tau }}\rangle$ on each site, and the numbers indicate the bond energies in unit of J. The long-range spin order we observed agrees with the conjecture based on physical insights in [7].

Zoom In Zoom Out Reset image size

We prefer the shorter, more descriptive name of 'dimer phase' adopted here because it indicates a solid (crystal) order of 'dimers', i.e. antiferromagnetically aligned spins along the bond direction, on a subset of the bonds. We propose to describe the long-range order using the vorticity or the winding number of the spin configuration around each hexagon

$$\begin{eqnarray}&&\nu =\displaystyle \frac{1}{3}\displaystyle \sum _{j=1}^{6}{\hat{n}}_{j}\cdot {\hat{n}}_{j+1},\end{eqnarray} \tag{ 15 }$$

where j labels the six sites of the hexagon. For example, hexagons marked by a dot in the center in figure 5(a) correspond to $\nu =1$ where all spins on the vertices point radially outwards (corresponding to the 'loop' in [7]). The rest of the hexagons, each marked by a cross at the center, have $\nu =-1/2$. It then becomes apparent that the hexagons marked by dots form a triangular lattice of spin vortices. And within one unit cell of the triangular lattice, the total vorticity is zero. Note that the state shown in figure 5(a) is energetically degenerate with a state where all the spins are flipped.

By embedding the quantum ${120}^{\circ }$ model into the more general tripod model, we are able to monitor the suppression of the dimer order and its eventual transition into the gapless SL phase (phase B) of the Kitaev model. The results are summarized in figure 6. We observe that the in-plane magnetization O₂ decreases continuously to zero as ${\theta }^{\prime }$ is reduced. Meanwhile, the ground state energy E steadily rises, indicating an increased degree of frustration as the Kitaev point is approached. One can measure the dimer order by introducing the energy difference $\delta E$ between the averages of two types of bonds: the 'happy' bonds (dimers) with antiparallel spins and the frustrated bonds where the two spins form an angle of 60°. Therefore, the dimensionless parameter

$$\begin{eqnarray}&&\eta =\delta E/E\end{eqnarray} \tag{ 16 }$$

can also serve as the order parameter for the dimer phase. As plotted in figure 6, η continuously drops to zero as ${\theta }^{\prime }$ is reduced from ${120}^{\circ }$ to 95°. Once inside the SL phase, both O₂ and η vanish, and the bond energies become approximately the same (see figure 5(b)). One can view the transition from the SL to the dimer phase as condensation of spin vortices. Equivalently, when ${\theta }^{\prime }$ is reduced, one can view the demise of the dimer order as the melting of the spin vortex lattice. Note the bond energy shown in figure 5 features small fluctuations and does not strictly obey C₆ rotation symmetry. In our tensor network calculations, no spatially symmetry is enforced on the tensors, and the expectation values of operators are computed approximately. The fluctuations are expected to decrease as the bond dimension is increased.

Zoom In Zoom Out Reset image size

The tripod model can be realized from the following Hubbard model at half filling in the Mott limit

$$\begin{eqnarray}&&{H}_{{\rm{hub}}}=-\displaystyle \sum _{i,\gamma ,\sigma {\sigma }^{\prime }}{t}_{\sigma {\sigma }^{\prime }}^{\gamma }{f}_{i,\sigma }^{\dagger }{f}_{i+{e}_{\gamma },{\sigma }^{\prime }}+U\displaystyle \sum _{i}{n}_{i,\uparrow }{n}_{i,\downarrow },\end{eqnarray} \tag{ 17 }$$

where ${f}_{i,\sigma }$ annihilates a fermion with spin σ at site i. The direction and spin dependent hopping matrix is related to ${T}^{\gamma }$ defined in equation (2) simply by

$$\begin{eqnarray}&&{t}_{\sigma {\sigma }^{\prime }}^{\gamma }=t{[1/2+{T}^{\gamma }(\theta )]}_{\sigma {\sigma }^{\prime }}.\end{eqnarray} \tag{ 18 }$$

Explicitly, they are given by

$$\begin{eqnarray}{t}_{\downarrow \downarrow } & = & \displaystyle \frac{t}{2}(1-{c}_{1}{c}_{2}),\;\;\;{t}_{\uparrow \uparrow }=\displaystyle \frac{t}{2}(1+{c}_{1}{c}_{2}),\\ {t}_{\uparrow \downarrow } & = & \displaystyle \frac{t}{2}({c}_{1}{s}_{2}-{\rm{i}}{s}_{1}),{t}_{\downarrow \uparrow }=\displaystyle \frac{t}{2}({c}_{1}{s}_{2}+{\rm{i}}{s}_{1}),\end{eqnarray}$$

where we have suppressed the superscript γ for brevity, and

$$\begin{eqnarray}&&{s}_{1}=\mathrm{sin}\theta ,{c}_{1}=\mathrm{cos}\theta ,{s}_{2}=\mathrm{sin}{\phi }_{\gamma },{c}_{2}=\mathrm{cos}{\phi }_{\gamma }.\end{eqnarray}$$

In the limit of $U\gg t$, using second-order perturbation theory, we obtain the following effective Hamiltonian for ${H}_{\mathrm{hub}}$

$$\begin{eqnarray}&&{H}_{\mathrm{eff}}=J\displaystyle \sum _{{\bf{r}},\gamma }{T}_{{\bf{r}}}^{\gamma }(\theta ){T}_{{\bf{r}}+{{\bf{e}}}_{\gamma }}^{\gamma }(\theta )-\displaystyle \frac{J}{4},\end{eqnarray} \tag{ 19 }$$

which is nothing but the tripod model $H(\theta )$, up to a constant term, with the superexchange $J=2{t}^{2}/U$. Note that the derivation of the effective compass Hamiltonian above does not depend on the details of the parameterization of ${T}^{\gamma }$ in terms of θ and ${\phi }_{\gamma }$. It follows that a large class of compass models, not limited to the tripod model proposed here, can be engineered on honeycomb lattice following the recipe above.

Duan et al previously showed that the Kitaev model can be realized using cold atoms on a hexagonal optical lattice with extra laser beams [36]. Generalization of their idea to the case of the tripod model (and other compass models) requires spin-dependent hopping ${t}_{\sigma {\sigma }^{\prime }}$ controlled by a non-Abelian gauge field or generalized spin–orbit coupling. Schemes to realize spin–orbit coupling was proposed in various approaches [37–40]. The realization of many have been demonstrated successfully in cold atoms experiments [41–44]. For example, spin-dependent optical lattices have been engineered using magnetic gradient modulation [44–46]. It seems possible, but challenging, to make ${t}_{\sigma {\sigma }^{\prime }}$ spatially dependent. Alternatively, the tripod model proposed here may be emulated using other artificial quantum systems such as superconducting quantum circuits [47].

The tripod model introduced in this paper encompasses three well known models of quantum magnetism: the Ising model, the Kitaev model and the 120° model. We established its (zero temperature) phase diagram using tensor network algorithms. This amounts to solving a continuum of frustrated spin models with spatially dependent exchange interactions. In particular, we found an extended SL phase around the Kitaev point, and a dimer phase for large values of angle ${\theta }^{\prime }$ including the quantum 120° model. The two quantum critical points obtained from tensor network states agree roughly with estimations from spin wave theory.

Our work only scratches the surface of the rich physics contained in the tripod model. Here we mention just a few open questions to be addressed in future work. First of all, it is desirable to develop a fieldtheoretical description of the continuous phase transition between the SL phase and the long-range ordered dimer phase, based on the intuitive picture of spin vortices introduced in section 5. Secondly, the tripod model, like other compass models, has very interesting emergent symmetry properties including intermediate symmetries midway between the global and local symmetries [1]. Consequently, the excitation spectrum is expected to contain zero modes and/or flat bands. It is therefore valuable to understand the excitation spectra of the long-range order phases by going beyond the ground state analysis here. Thirdly, the finite temperature properties of the tripod model deserve a separate study. The classical limit of the tripod model is known to be highly nontrivial. The effects of thermal fluctuations and the 'order by disorder' mechanism have been investigated in [34] for the classical 120° model. Finally, we have only focused on the case of ${J}_{\gamma }=J$ here. From the Kitaev model, we know that a gapped SL phase (phase A) takes over when the asymmetry in J_γ grows large. Thus one expects that further generalization of the tripod model to general values of J_γ may uncover new interesting phases.

To conclude, we hope our results can stimulate further application of tensor network algorithms to frustrated spin models as well as spin–orbital models describing transition metal oxides. We also hope our introduction of the tripod model can inspire alternative proposals to extend the Kitaev model or realize compass models in artificial quantum systems such as cold atoms on optical lattices or superconducting circuits.

We thank Jiyao Chen, Arun Paramekanti, Zhiyuan Xie, and Congjun Wu for helpful discussions. This work is supported by AFOSR Grant No. FA9550-16-1-0006 (HZ, BL, EZ, and WVL), NSF Grant No. PHY-1205504 (HZ and EZ), and jointly by US ARO Grant No. W911NF-11-1-0230, the Pittsburgh Foundation and its Charles E Kaufman Foundation, and the Overseas Scholar Collaborative Program of NSF of China No. 11429402 sponsored by Peking University (WVL). WVL is grateful for the hospitality of KITP UCSB where part of this manuscript was written with support from NSF PHY11-25915. Publication of this article is funded by the George Mason University Libraries Open Access Publishing Fund, and the University Library System, University of Pittsburgh.

The finite PEPS algorithm is a powerful numerical approach for two-dimensional quantum spin systems [25–27]. For the tripod model, we construct the usual PEPS wave function starting from six random rank-four tensors with virtual bond dimension D = 3. The tensors are then optimized through recursive imaginary-time evolution with time step $\tau =0.01$ on all the links. Once the wave function is converged, we calculate the ground state energy and the expectation value of the combined order parameter $O=\sqrt{{O}_{1}^{2}+{O}_{2}^{2}}$. The results are shown in figure 7 for three typical values of ${\theta }^{\prime }$ (corresponding to the three different phases found in the thermodynamic limit). For the small system size considered here (six sites with periodic boundary conditions), O vanishes as the wave function converges to the ground state. Nonetheless, a noticeable peak of O during the time evolution can serve as the indication for spin order. For all the three cases, the energies converge to the ED result up to a relative error of 10⁻³ (the upper panel of figure 7). The peaks of O in the lower panel of figure 7 for the cases ${\theta }^{\prime }=70^\circ$ and ${\theta }^{\prime }=120^\circ$ point to the Néel order phase and the dimer phase respectively, while the monotonic decay of O for the case ${\theta }^{\prime }=90^\circ$ suggests a SL state.

Zoom In Zoom Out Reset image size

Within the many variants of tensor network algorithms, a typical way to find the phase diagram of quantum spin models is the infinite PEPS (iPEPS) method [28, 29]. The iPEPS ansatz on the honeycomb lattice usually proceeds by mapping the lattice to a square lattice and evaluating the effective environment by contraction schemes such as infinite matrix product states [28] or corner transfer matrices [48, 49]. For instance, the phases of Kitaev–Heisenberg model [50] and the SU(4) symmetric Kugel–Khomskii model [51] have been studied via the iPEPS ansatz with a 2 × 2 or 4 × 4 unit cell. However, the contraction scheme for a six-site (hexagonal) unit cell on the honeycomb lattice is tedious and expensive, especially for the corner transfer matrices scheme.

For this reason, we adopt the simple tensor update scheme and evaluate the contraction using the HOTRG method as explained in the main text. The simple update [32] generalizes the TEBD [33] technique to two-dimensional quantum systems by introducing the bond vectors to represent the mean-field environment for local tensors. We set the imaginary time step $\tau =0.01$ and the number of iterations is generally around 10⁵ (smaller time step does not improve the numerical result significantly). The accuracy of the HOTRG method is controlled by the virtual bond dimension D. By systematically increasing D, the quantum entanglement between neighboring sites is better taken into account, yielding a more accurate ground state. For example, figure 8 shows that the order parameter O₂ vanishes when D is increased to 8 in the region ${\theta }^{\prime }\lt 94^\circ$, suggesting a SL ground state. One notices that the variations of the ground state energy with ${\theta }^{\prime }$ within the SL phase is larger than those in the long-range ordered phase, especially for smaller D values. This is due to the strong quantum fluctuations intrinsic to the SL. Cross-checking the HOTRG calculations here to those using Second Renormalization Group [52] which takes into account the entanglement between the system and the environment deserves a future study.

Zoom In Zoom Out Reset image size

For the two long-range ordered phases, the $1/S$ corrections to the energy are obtained by diagonalizing the matrix equation (14) with different form factors ${\beta }_{1},{\beta }_{2}$, and ${\beta }_{3}$. For the Néel ordered phase, they are given by

$$\begin{eqnarray}{\beta }_{1} & = & \displaystyle \sum _{j=1}^{3}{({b}_{j}+{\rm{i}}{a}_{j})}^{2}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\hat{e}}_{j}},\\ {\beta }_{2} & = & \displaystyle \sum _{j=1}^{3}{({b}_{j}+{\rm{i}}{a}_{j})}^{2}{{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\hat{e}}_{j}},\\ {\beta }_{3} & = & -\displaystyle \sum _{j=1}^{3}({b}_{j}^{2}+{a}_{j}^{2}){{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\hat{e}}_{j}}.\end{eqnarray}$$

And for the dimer phase

$$\begin{eqnarray}{\beta }_{1} & = & \displaystyle \sum _{j=1}^{3}{({a}_{j}+{\rm{i}}{b}_{j})}^{2}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\hat{e}}_{j}},\\ {\beta }_{2} & = & \displaystyle \sum _{j=1}^{3}{({a}_{j}+{\rm{i}}{b}_{j})}^{2}{e}^{-i{\bf{k}}\cdot {\hat{e}}_{j}},\\ {\beta }_{3} & = & -\displaystyle \sum _{j=1}^{3}({a}_{j}^{2}+{b}_{j}^{2}){{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\hat{e}}_{j}}.\end{eqnarray}$$

Here, a_j and b_j are related to the parameter $\theta ,\varphi$, and α defined in the main text through

$$\begin{eqnarray}{a}_{1} & = & \mathrm{cos}\theta \mathrm{cos}\varphi \mathrm{cos}\alpha -\mathrm{sin}\theta \mathrm{sin}\varphi ,\\ {a}_{2} & = & \mathrm{cos}\theta \mathrm{cos}\varphi \mathrm{cos}\left(\alpha -\displaystyle \frac{2\pi }{3}\right)-\mathrm{sin}\theta \mathrm{sin}\varphi ,\\ {a}_{3} & = & \mathrm{cos}\theta \mathrm{cos}\varphi \mathrm{cos}\left(\alpha -\displaystyle \frac{4\pi }{3}\right)-\mathrm{sin}\theta \mathrm{sin}\varphi ,\\ {b}_{1} & = & -\mathrm{cos}\theta \mathrm{sin}\alpha ,\\ {b}_{2} & = & -\mathrm{cos}\theta \mathrm{sin}\left(\alpha -\displaystyle \frac{2\pi }{3}\right),\\ {b}_{3} & = & -\mathrm{cos}\theta \mathrm{sin}\left(\alpha -\displaystyle \frac{4\pi }{3}\right).\end{eqnarray}$$