A continuum of compass spin models on the honeycomb lattice

Quantum spin models with spatially dependent interactions, known as compass models, play an important role in the study of frustrated quantum magnetism. One example is the Kitaev model on the honeycomb lattice with spin-liquid ground states and anyonic excitations. Another example is the geometrically frustrated quantum $120^\circ$ model on the same lattice whose ground state has not been unambiguously established. To generalize the Kitaev model beyond the exactly solvable limit and connect it with other compass models, we propose a new model, dubbed"the tripod model", which contains a continuum of compass-type models. It smoothly interpolates the Ising model, the Kitaev model, and the quantum $120^\circ$ model by tuning a single parameter $\theta'$, the angle between the three legs of a tripod in the spin space. Hence it not only unifies three paradigmatic spin models, but also enables the study of their quantum phase transitions. We obtain the phase diagram of the tripod model numerically by tensor networks in the thermodynamic limit. We show that the ground state of the quantum $120^\circ$ model has long-range dimer order. Moreover, we find an extended spin-disordered (spin-liquid) phase between the dimer phase and an antiferromagnetic phase. The unification and solution of a continuum of frustrated spin models as outline here may be useful to exploring new domains of other quantum spin or orbital models.


I. INTRODUCTION
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][3][4][5][6][7] as well as interacting systems with spinorbit 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 Ref. [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 σ x i σ x j , σ y i σ y j , and σ z i σ z j 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 * ezhao2@gmu.edu 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 spin liquid 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][15][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 σ x , σ y , and σ z for the three bond directions are replaced by three (pseudo)spin 1/2 operators T 1 , T 2 , and T 3 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. orbitalonly 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 θ and a simple, intuitive picture for the three (pseudo)spin operators: they form a tripod in spin space as shown in Fig. 1. Analogous to tuning the tripod to raise or low a mounted camera in photography, dialing the angle θ 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 spin liquids.
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][20][21][22]. The results are summarized in Fig. 1 and Fig. 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 spin liquid phase exists in an extended region in our phase diagram [25]. The numerical results of TRG are further confirmed and crosschecked with Projected Entangled Pair States (PEPS) calculations [26,27] for finite systems, exact diagonalization, and spin wave analysis. We discuss the qualitative features of the quantum phase transitions between the spin liquid 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.

II. THE TRIPOD MODEL
We generalize the Kitaev and the 120 • model to the following continuum compass model defined on the twodimensional honeycomb lattice, where J > 0, and for each lattice site r the spin 1/2 operators are defined as with the Pauli matrices τ x , τ y , τ z . Each site r is coupled to its neighbors r + e γ , where e γ , γ = 1, 2, 3, denotes the three bond vectors of the honeycomb lattice. Geometrically, the three T γ form a tripod in the spin space as shown in Fig. 1 (a) The nearest neighbor spin exchange along bond direction eγ, γ = 1, 2, 3, is defined through the spin 1/2 operator T γ , represented by an arrow in spin space spanned by τ x , τ y , τ z . The three T γ can be thought as the three legs of a tripod, being tilted out of the xz plane by angle θ, and forming an angle θ with each other. When projected onto the xz plane, T γ is along the eγ direction. (b) The schematic phase diagram of the tripod model. As the tripod is opened by increasing θ , the model starts as the Ising model at θ = 0, becomes the Kitaev model at θ = 90 • , and then the quantum 120 • model at θ = 120 • . Three phases are identified: a Néel ordered antiferromagnet, a spin liquid (SL), and a long-range ordered dimer phase.
angle θ and, when projected onto the xz plane, are orientated along the corresponding bond direction e γ , i.e. at azimuthal angle φ γ = 0, 2π/3, 4π/3 respectively. While T γ is most naturally defined through the tilting angle θ, it is much more convenient to introduce another angle, θ , to discuss the various limits of H(θ). θ is the angle between T 1 and T 2 , i.e., the two adjacent legs of the tripod. And it is related to θ by trigonometry We will take θ as the only tuning parameter in the tripod model. Three special limits of this model can now be identified. First of all, when θ = 0, T γ all collapse to τ y . The tripod is closed, and H(θ) is nothing but the Ising model, Note that we choose τ y as the vertical axis in spin space instead of the usual convention of τ z so that H(θ) reduces exactly to the 120 • model defined in our earlier work Ref. [6]. Secondly, when θ = 90 • , H(θ) reduces to the Kitaev model since the operators T γ are now perpendicular to each other in the spin space. We can simply identify them as σ x , σ y and σ z (apart from a factor 1/2) in a rotated coordinate system as illustrated in Fig. 1(b). Thirdly, for θ = 120 • , H(θ) 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 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 By solving H(θ) using tensor network algorithms, we compute the average spin τ x,y,z r in the ground state. The main results are summarized in the schematic phase diagram in Fig. 1(b). Fig. 2 shows the variation of the two order parameters introduced above as θ is changed. The region at small θ corresponds to the familiar Néel order which is characteristic of the classical Ising model and illustrated in the left inset of Fig. 2. Despite the increased quantum fluctuations as θ is increased, the Néel ordered phase persists up to θ ∼ 87 • . At the opposite end of large θ , 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 [ Fig. 5(a)]. The triangular lattice and its enlarged unit cell becomes transparent if we introduce a topological charge [red dot in Fig. 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 Fig. 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 • model on honeycomb lattice. This point will be further discussed in Section V.
Sandwiched between the Néel ordered phase and the dimer phase, a quantum spin liquid phase is stabilized for θ ∈ [87 • , 94 • ]. The conclusion is mainly based on the observation from Fig. 2 that the order parameters O 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 θ = 90 • . The order parameters as functions of θ in Fig. 2 also suggest that the two quantum phase transitions in the tripod model may be qualitatively different. The gradual drop of O 2 at θ > 90 • indicates a continuous phase transition between the dimer phase and the spin liquid phase. In contrast, the drop of the order parameter O 1 at θ < 90 • 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 III.

III. TENSOR NETWORK ALGORITHMS
Recent developments of entanglement-based tensor network algorithms provide a novel, accurate approach to strongly correlated electron systems [26][27][28][29][30]. Particularly, they have been successfully applied to frustrated quantum magnets [19][20][21][22] and the t − J model [31,32] 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 [26][27][28] to solve H(θ) for a six-site system with periodic boundary conditions. The ground state energies obtained coincide with those from exact diagonalization. This suggests that PEPS is intrinsically superior compared to mean field theories when applied to frustrated spin Hamiltonians such as H(θ). The order parameters decay to zero as the ground state is approached for a finite system. Nonetheless, their decay behaviors are quite disparate for θ values in the Ising, Kitaev, and 120 • 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 Ref. [33] which generalizes the time-evolving block decimation (TEBD) [34] technique to two dimensions. For a n-site unit cell, e.g. a six-site unit cell shown in Fig. 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 |Ψ is characterize by n tensors T j , j ∈ [1, n]. We then evaluate the expectation value of operator O, O = Ψ|O|Ψ / Ψ|Ψ which involves the (infinite) product of tensors T j , using a real space coarse graining procedure known as higher-order tensor renormalization group (HOTRG) [24] schematically shown in Fig. 3. We outline the main steps here. At the i-th step of HOTRG, a local tensor, say T i 1 , is regrouped with its three nearest neighbor tensors (T 2 , T 4 , T 6 ) to form a new tensorT i+1 1 . More generally, for odd or even sites, T i+1 e = s.l.
where the summation is over the shared legs [abbreviated as s.l., the solid lines in Fig. 3 where the three projection tensors U x,y,z , shown in Fig. 3(b), are obtained as follows. Take U x as an example. First,T j is reshaped into matrix T j with the row corresponding to the leg along the x-direction. Then, a matrix U is obtained by singular value decomposition (SVD) Finally, U x is obtained by truncating 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 i but represent a larger system, see Fig. 3(c). This constitutes a single RG step. By iterating the RG steps many times, the converged result of O well approximates the expectation value in the thermodynamic limit.
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 θ < 90 • . 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 θ values in the ground state calculation. In practice, one can safely use the two-site unit cell for θ < 90 • since it is significantly cheaper. The phase diagram and the spin configurations in the ordered phases shown in Fig. 2 are obtained by using the two-site unit cell for θ ≤ 90 • and the six-site unit cell for θ > 90 • .

IV. SPIN WAVE ANALYSIS
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 r = ±T r for all sites on the A (B) sublattice. Then the tripod Hamiltonian acquires a suggestive form Here we have promoted the spin 1/2 operator τ /2 to general spin operator 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 γ r (θ) = S γ r+eγ (θ), i.e., the projection of 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 Ref. [1]. The special case of the classical 120 • model on honeycomb lattice was previously discussed in Ref. [7,35]. We will confine our spin wave analysis to the simple case of spatially homogeneous spin configurations S r = S 0 as done in Ref. [35]. The direction of S 0 is characterized by its polar angle ϕ measured from τ y and its azimuthal angle α of S 0 measured from τ z in τ x -τ z plane. The corresponding classical ground state energy per unit cell is (2 sin 2 θ cos 2 ϕ + cos 2 θ sin 2 ϕ).
It is interesting to note that, coincidently, at the Kitaev point, θ = 90 • which corresponds to θ = cos −1 2/3, E 0 is completely flat and does not depend on ϕ. For θ < 90 • , E 0 is minimized when ϕ = 0 or π, corresponding to the two degenerate states with spin up or down in the Néel ordered phase. In contrast, for θ > 90 • , E 0 is minimized when ϕ = π/2, i.e., S 0 lies within the τ x -τ 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 [36] to H(θ) 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 where ∆ = −E 0 /S 2 J, N is the number of sites within the A sublattice, the k summation is over the first Brillouin zone of the A sublattice, and ω λ (k) describe the spin wave dispersion for branch λ = 1, 2 and they are given by the eigenvalues of the matrix The expressions for β 1 , β 2 and β 3 are lengthy and tabulated in Appendix B. In what follows, we will discuss the energy E(ϕ, α) = E 0 + E 1 separately for the two distinct cases: θ < 90 • and θ > 90 • . The results for E(ϕ) from the spin wave analysis are plotted in the upper panel of Fig. 4 for several values of θ corresponding to the Néel ordered phase. One notices that the fluctuations do not change qualitatively the mean field ground state. E(ϕ) reaches minima still at ϕ = 0 or π for small θ . However, as θ is increased, the energy E(ϕ) becomes flatter. Eventually, as θ = θ c ∼ 75.0 • (the top curve of Fig. 4), the energies for ϕ = π/4, 3π/4 with proper choice of α become degenerate with those for ϕ = 0, π. This signals the destabilization of the Néel order by quantum fluctuations. This occurs around θ c , before the Kitaev point is approached. Note that in Fig. 4, only the region ϕ ∈ [0, π/4] ∪ [3π/4, π] is shown. Outside this region (and also for θ > θ c ), the lowest order spin wave theory based on the Ising-like antiferromagnetic order becomes ill defined.
For θ > 90 • , the classical ground state is continuously degenerate with ϕ = π/2 but arbitrary α ∈ [0, 2π]. 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 θ = 120 • , i.e. the quantum 120 • model, has been discussed before in Ref. [6,7,35]. As shown in the lower panel of Fig. 4, the same physical picture continues to hold for the tripod model for θ ≤ 120 • : the energy E is minimized at α n = nπ/3 with integer n. However, E(α) becomes increasingly flat as θ is decreased. At the critical point θ = θ d ∼ 94.6 • , additional minima of E appear at α = α n +π/6. This indicates that the long-range order gets destroyed and replaced by a new phase for θ ≤ θ d .
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 spin liquid phase, θ c and θ d . 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 θ > 90 • .

V. THE DIMER PHASE
Previous theoretical studies of the quantum 120 • model on the honeycomb lattice gave conflicted results. The spin wave analysis of Ref. [6] assumed a homogenous ground state S r = S 0 and found quantum fluctuations prefer α n = nπ/3. This led the authors to suggest that the ground state may be a simple ferromagnet of S with any choice of α n (i.e. antiferromagnetic order in terms of the original spin T or τ ). Ref. [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 α 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 • model has a ground state with the six-site plaquette order [7]. However, no evidence of long-range order was found in exact diagonalization (ED) studies where the spin correlation functions were computed for finite size clusters with periodic boundary conditions [35]. Instead, the ED results supported a trial wave function similar in spirit to the shortrange 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 • 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 • model is the dimer phase illustrated in Fig. 5(a) where the arrows denote the directionn of spin average τ 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 Ref. [7].
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, where j labels the six sites of the hexagon. For example, hexagons marked by a dot in the center in Fig. 5(a) correspond to ν = 1 where all spins on the vertices point radially outwards (corresponding to the "loop" in Ref. [7]). The rest of the hexagons, each marked by a cross at the center, have ν = −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 Fig. 5(a) is energetically degenerate with a state where all the spins are flipped. By embedding the quantum 120 • 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 spin liquid phase (phase B) of the Kitaev model. The results are summarized in Fig. 6. We observe that the in-plane magnetization O 2 decreases continuously to zero as θ 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 δ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 can also serve as the order parameter for the dimer phase. As plotted in Fig. 6, η continuously drops to zero as θ is reduced from 120 • to 95 • . Once inside the spin liquid phase, both O 2 and η vanish, and the bond energies become approximately the same [see Fig. 5(b)]. One can view the transition from the spin liquid to the dimer phase as condensation of spin vortices. Equivalently, when θ 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 Fig. 5 features small fluctuations and does not strictly obey C 6 rotation symmetry. In our tensor network calculations, no spatially symmetry is enforced on the tensors, and the

VI. POTENTIAL REALIZATION
The tripod model can be realized from the following Hubbard model at half filling in the Mott limit, where f i,σ annihilates a fermion with spin σ at site i. The direction and spin dependent hopping matrix is related to T γ defined in Eq. (2) simply by Explicitly, they are given by where we have suppressed the superscript γ for brevity, and s 1 = sin θ, c 1 = cos θ, s 2 = sin φ γ , c 2 = cos φ γ .
In the limit of U t, using second-order perturbation theory, we obtain the following effective Hamiltonian for H hub which is nothing but the tripod model H(θ), up to a constant term, with the superexchange J = 2t 2 /U . Note that the derivation of the effective compass Hamiltonian above does not depend on the details of the parameterization of T γ in terms of θ and φ γ . 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 [37]. Generalization of their idea to the case of the tripod model (and other compass models) requires spin-dependent hopping t σσ controlled by a non-Abelian gauge field or generalized spin-orbit coupling. Schemes to realize spin-orbit coupling was proposed in various approaches [38][39][40][41]. The realization of many have been demonstrated successfully in cold atoms experiments [42][43][44][45]. For example, spindependent optical lattices have been engineered using magnetic gradient modulation [45][46][47]. It seems possible, but challenging, to make t σσ spatially dependent. Alternatively, the tripod model proposed here may be emulated using other artificial quantum systems such as superconducting quantum circuits [48].

VII. OUTLOOK
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 spin liquid phase around the Kitaev point, and a dimer phase for large values of angle θ 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 field-theoretical description of the continuous phase transition between the spin liquid phase and the long-range ordered dimer phase, based on the intuitive picture of spin vortices introduced in Section V. 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 Ref. [35] for the classical 120 • model. Finally, we have only focused on the case of J γ = J here. From the Kitaev model, we know that a gapped spin liquid 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. The finite PEPS algorithm is a powerful numerical approach for two-dimensional quantum spin systems [26][27][28]. 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 τ = 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 The results are shown in Fig. 7 for three typical values of θ (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 exact diagonalization result up to a relative error of 10 −3 (the upper panel of Fig. 7). The peaks of O in the lower panel of Fig. 7 for the cases 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 [29,30]. 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 [29] or corner transfer matrices [49,50]. For instance, the phases of Kitaev-Heisenberg model [51] and the SU(4) symmetric Kugel-Khomskii model [52] 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 higherorder tensor renormalization group (HOTRG) method as explained in the main text. The simple update [33] generalizes the time-evolving block decimation [34] 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 τ = 0.01 and the number of iterations is generally around 10 5 (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, Fig. 8 shows that the order parameter O 2 vanishes when D is increased to 8 in the region θ < 94 • , suggesting a spin liquid ground state. One notices that the variations of the ground state energy with θ within the spin liquid 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 spin liquid. Cross-checking the HOTRG calculations here to those using Second Renormalization Group [53] which takes into account the entanglement between the system and the environment deserves a future study. β 3 . For the Néel ordered phase, they are given by And for the dimer phase, (a j + ib j ) 2 e −ik·êj , Here, a j and b j are related to the parameter θ, ϕ, and α defined in the main text through a 1 = cos θ cos ϕ cos α − sin θ sin ϕ, a 2 = cos θ cos ϕ cos(α −