Exact chiral spin liquids and mean-field perturbations of gamma matrix models on the ruby lattice

The study of zero-temperature properties of 'frustrated' magnetic systems has been reinvigorated in recent years with the development of powerful numerical methods [1–3], the discovery of new classes of exactly solvable models [4, 5] and the realization that intriguing topological ground states could occur [6, 7]. In particular, there is a large class of exactly solvable quantum spin models known as Kitaev models [5] which have revitalized the study of exactly solvable spin-liquid systems. The appeal of this class of models lies in the relative ease with which Kitaev's original version can be generalized, spawning many variants [8–19], and its possession of many nontrivial properties [20–29] even in the presence of disorder [30–32]. Moreover, the exact solution of Kitaev models does not require any more effort than solving the problem of noninteracting particles moving in a static background magnetic field. Nonetheless, the exact eigenstates of all Kitaev models are nontrivial entangled many-body wavefunctions [22, 23]. More specifically, their respective ground states are examples of quantum spin liquids, which are insulating quantum states of matter exhibiting no conventional long-range magnetic order at zero temperature. Although no physical example of these models has been established in nature, there have been several experimental proposals for realizing them [33–39]. Nevertheless, the importance of this entire class of models lies in providing a strong case—by proof of principle—for the existence of exotic emergent phenomena in many-body quantum phases of matter such as quantum spin liquids, as well as their utility as model systems for the study of nontrivial, nonperturbative emergent quantum phenomena in general.

In this paper, we study a new exactly solvable Kitaev model and examine its response to the inclusion of interactions that destroy the exact solvability. This will motivate us to analyze the order that is favored by these new interactions and their effects on the exactly solvable ground state. The model is a gamma matrix model (GMM) extension [11, 12, 16] on the two-dimensional ruby lattice. Previous studies of the ruby lattice have yielded a topological insulator [40], fractional quantum anomalous Hall states [41] and topological anyons [42]. In this paper, we show that a GMM on the ruby lattice realizes a quantum spin liquid ground state with an unusually rich phase diagram. The model is then generalized by the inclusion of Ising-like interactions that spoil the exact solvability of the model. We explore the interplay between these interactions and the many ground state spin-liquid phases.

This paper is organized as follows. In section 2, we introduce the Hamiltonian and discuss its solution by mapping it onto a system of noninteracting Majoranas moving in an emergent $\mathbb {Z}_2$ gauge field. In section 3, we discuss the various ground state phases that the model realizes and present a phase diagram. In section 4, we examine the Ising order in one of the gamma matrix operators which may be induced by including Ising-like interactions. In section 5, we treat these extra interactions within the mean-field approximation and discuss their effects on the phase diagram. We then present the conclusions in section 6.

We consider spin-3/2 moments located at the sites of a ruby lattice, shown in figure 1(a). Alternatively, one can think of each site as being occupied by a single spin-1/2 which has an additional orbital degree of freedom. At each site, we introduce five 4 × 4 Hermitian gamma matrices which satisfy a Clifford algebra locally {Γ^a_i,Γ^b_i} = 2δ_ab, with $a,b=1,\dots ,5$ [11, 12, 16]. These matrices also have a representation in terms of bilinears of SU(2) spin-3/2 operators on each site i as follows:

$$\begin{eqnarray} \Gamma_i^1 &=& \frac{1}{\sqrt{3}} \{ S_i^y , S_i^z \}, \quad \Gamma_i^2 = \frac{1}{\sqrt{3}} \{ S_i^x,S_i^z \}, \nonumber \\ \Gamma_i^3 &=& \frac{1}{\sqrt{3}} \{ S_i^x , S_i^y \}, \quad \Gamma_i^4 = \frac{1}{\sqrt{3}} [(S_i^x)^2 - (S_i^y)^2] , \\ \Gamma_i^5 &=& (S_i^z)^2 - \textstyle\frac{5}{4}.\nonumber \end{eqnarray} \tag{ 1 }$$

The gamma matrices may also be expressed as tensor products of Pauli spin-1/2 matrices {σ^μ,τ^ν} on the site i,

$$\begin{eqnarray} \Gamma_i^1 &=& \sigma_i^z\otimes \tau_i^y, \quad \Gamma_i^2 = \sigma_i^z\otimes\tau_i^x, \nonumber \\ \Gamma_i^3 &=& \sigma_i^y\otimes\tau_i^0, \quad \Gamma_i^4 = \sigma_i^x\otimes\tau_i^0, \\ \Gamma_i^5 &=& \sigma_i^z\otimes\tau_i^z,\nonumber \end{eqnarray} \tag{ 2 }$$

where σ⁰_i = τ⁰_i = 1, the 2 × 2 identity matrix. In the spin–orbit interpretation (2), σ acts on the real spin degree of freedom and τ on the orbital degree of freedom. Thus, one can view the GMM we study below as applying to both a spin-3/2 model and a two-orbital spin-1/2 model. In either case, there are four physical states per lattice site.

Zoom In Zoom Out Reset image size

The Hamiltonian we study is given by

$$\begin{eqnarray} &&\fl \mathcal{H}= \mathcal{H}_0 + \mathcal{H}_5, \nonumber \\ &&\fl \mathcal{H}_0 = J_{\triangle} \sum_{\langle ij \rangle \in \triangle} \Gamma^1_i \Gamma^2_j + J'_{\triangle} \sum_{\langle ij \rangle \in \triangle} \Gamma^{15}_i \Gamma^{25}_j + J_{\unicode{x2394}} \sum_{\langle ij \rangle \in \unicode{x2394}} \Gamma^3_i \Gamma^4_j + J'_{\unicode{x2394}} \sum_{\langle ij \rangle \in \unicode{x2394}} \Gamma^{35}_i \Gamma^{45}_j, \\ &&\fl \mathcal{H}_5=- J_5 \sum_i \Gamma^5_i,\nonumber \end{eqnarray} \tag{ 3 }$$

where we define the operators Γ^ab_i = [Γ^a_i,Γ^b_i]/(2i) in terms of the Γ^a_i given in either equation (1) or equation (2). The couplings $J_{\triangle },J'_{\triangle },J_{\unicode{x2394} }$ and $J'_{\unicode{x2394} }$ are taken to be positive and are link dependent, as shown in figure 1(b). For generic couplings, the Hamiltonian has translational and sixfold rotational lattice symmetry. With the interpretation in terms of spin-3/2 moments, it has global Ising spin symmetry under 180° rotations about the z-axis, and time-reversal symmetry (TRS), although TRS will be spontaneously broken in the ground state as we will be show below.

While the model (3) has unusual spin symmetries in terms of the underling spin and orbital degrees of freedom, its structure allows an exact solution. A key ingredient for the exact solvability of the model is the existence of conserved operators on every plaquette of the lattice which commute with the Hamiltonian [5]. These operators are the three ($\hat {W}_{\triangle }$), four ($\hat {W}_\square$) and six point ($\hat {W}_{\unicode{x2394} }$) operators defined on the triangles, squares and hexagons of the lattice (as their subscripts suggest). Explicitly, they are given by $\hat {W}_{\triangle } = \Gamma ^{12}_i \Gamma ^{12}_j \Gamma ^{12}_k,$ $\hat {W}_{\unicode{x2394} } = \Gamma ^{34}_i \Gamma ^{34}_j \Gamma ^{34}_k \Gamma ^{34}_l \Gamma ^{34}_m \Gamma ^{34}_n$ and $\hat {W}_{\square } = \Gamma ^{23}_i \Gamma ^{14}_j \Gamma ^{23}_k \Gamma ^{14}_l,$ where the sites are labeled counter-clockwise and for $\hat {W}_{\square }$ the first link 〈ij〉 lies on a triangle. The model is then solved by introducing a Majorana representation of the Γ-matrices using six flavors of Majorana operators, {ξ¹_i,ξ²_i,ξ³_i,ξ⁴_i,c_i,d_i}, at each site i:

$$\begin{equation} \Gamma_{i}^{a} = \mathrm{i}\,\xi_{i}^{a}c_{i},\quad \Gamma_{i}^{5} = \mathrm{i}c_{i}d_{i},\quad \Gamma_{i}^{a5} = \mathrm{i}\,\xi_{i}^{a}d_{i}, \end{equation} \tag{ 4 }$$

where the Majorana fermions have the important property that (ξ^a_i)^† = ξ^a_i, c^†_i = c_i and d^†_i = d_i. The Hamiltonian written in terms of these operators then reduces to that of two species of Majorana fermions, c and d, which are noninteracting and move in a background $\mathbb {Z}_2$ gauge field u_ij,

$$\begin{eqnarray} &&\fl \tilde{\mathcal{H}} = J_{\triangle} \sum_{\langle ij \rangle\in\triangle} \mathrm{i} u_{ij} c_i c_j + J'_{\triangle} \sum_{\langle ij \rangle\in\triangle} \mathrm{i} u_{ij} d_i d_j + J_{\unicode{x2394}} \sum_{\langle ij \rangle\in\unicode{x2394}} \mathrm{i} u_{ij} c_i c_j + J'_{\unicode{x2394}} \sum_{\langle ij \rangle\in\unicode{x2394}} \mathrm{i} u_{ij} d_i d_j - J_5 \sum_i \mathrm{i} c_i d_i.\nonumber\\ \end{eqnarray} \tag{ 5 }$$

Here the fields u_ij are defined by u_ij = −i ξ¹_iξ²_j if $ij \;\in \triangle$ and u_ij = −i ξ³_iξ⁴_j if $ij\in \unicode{x2394}$. The u_ij are, in general, quantum fields with eigenvalues ±1 since u²_ij = 1, hence their identification with a $\mathbb {Z}_2$ gauge theory. Moreover, we can simultaneously diagonalize ${\tilde {\mathcal {H}}}$ and {u_ij} since they commute. We emphasize that the interacting Hamiltonian (3) with the representation (4) is reduced to the effectively noninteracting Hamiltonian (5) because the u_ij behave as constants in each flux sector of the $\hat {W}_{\triangle }$, $\hat {W}_\square$ and $\hat {W}_{\unicode{x2394} }$.

However, the full Hilbert space spanned by the Majorana fermions is overcomplete, and a constraint must be enforced to ensure that the Clifford algebra of Γ operators is satisfied. This constraint is expressed by the operator equation D_i = − Γ¹_iΓ²_iΓ³_iΓ⁴_iΓ⁵_i = − i ξ¹_iξ²_iξ³_iξ⁴_ic_id_i = 1 and is enforced by the projector $P=\prod _{i} [\frac {1+D_{i}}{2} ]$, where the product is over all sites in the lattice. The original Hamiltonian is obtained through ${\cal H}=P \tilde{{\cal H}}P$.

The Z₂ gauge fields define fluxes ϕ_p via $\exp (\mathrm {i}\phi _p)\equiv \prod _{jk\in p} \mathrm {i}u_{jk}$, where jk is taken counterclockwise on each elementary plaquette p. These fluxes are then related to the conserved $\hat {W}$s by $W_{p} \propto \prod _{ij\in p}u_{ij}$, where ij is also taken counterclockwise. Since the eigenstates of the Hamiltonian are also eigenstates of the fluxes, once the fluxes have been specified, the u_ij = ± 1s are uniquely specified up to $\mathbb {Z}_2$ gauge transformations, and the Hamiltonian ${\skew3\tilde{\mathcal {H}}}$ is then diagonalized and projected to yield to eigenstates of $\mathcal {H}$. The ground state is then determined by a many-body Majorana wavefunction minimizing the total energy. Due to translational symmetry, ${\skew3\tilde{\mathcal{H}}}$ describes a band structure for the c and d fermions. It can be shown that the minimal energy configuration is the one where all the negative 'eigenstates' of the effective band Hamiltonian (5) are occupied [5, 11, 12].

One finds that ϕ_p = ± π/2 in the triangular plaquettes, and ϕ_p = 0,π in the hexagonal and square plaquettes. Under TRS, W_p → ± W_p where −(+) is for triangle (hexagon and square) plaquettes; it follows that ϕ_p → − ϕ_p for triangle plaquettes, while ϕ_p remains unchanged for hexagon and square plaquettes. Consequently, a ground state with a certain flux pattern {ϕ_p} spontaneously breaks TRS. The ground state energy of a flux configuration must be degenerate with the flux pattern obtained from {ϕ_p} by changing ϕ_p → − ϕ_p on all triangular plaquettes [5, 11].

To determine the physical ground state, the flux configuration that minimizes the ground state energy needs to be determined. This was accomplished by first numerically determining the ground state energy for each flux configuration for the symmetric couplings $J_{\triangle }=J'_{\triangle }=J_{\unicode{x2394} }=J'_{\unicode{x2394} }=J$ and J₅ = 0. The ground state flux with the least energy is depicted by figure 1(b). Generally, the ground state flux will be a function of J₅, but for simplicity we will not consider additional flux configurations in this work because an additional term (depending on $\hat {W}_{\triangle },\hat {W}_{\square },\hat {W}_{\unicode{x2394} }$) can always be included to favor a particular configuration without destroying the exact solvability [20, 28]. By general arguments regarding the Majorana representation of the spin operators in the 3/2 representation [5, 25], the spin–spin correlation functions of the ground state are identically 0 beyond nearest-neighbor sites. Hence, the ground state is also a quantum spin liquid in these observables.

With the gauge sector of $\mathcal {H}$ specified, we can solve our system by mapping the Majorana fermions to complex fermions a_i, defined by c_i = a_i + a^†_i and d_i = −i(a_i − a^†_i). This puts our Hamiltonian into the form

$$\begin{equation} \fl {\tilde {\mathcal {H}}} = \sum_{\langle ij \rangle} \{ \mathrm{i} (u^c_{ij} + u^d_{ij}) a^\dagger_i a_j + \textstyle\frac{1}{2} i (u^c_{ij} - u^d_{ij}) (a^\dagger_i a^\dagger_j + a_i a_j) \}- 2 \hspace{.4mm} J_5 \sum_i (a^\dagger_i a_i - \textstyle\frac{1}{2}), \end{equation} \tag{ 6 }$$

where we have defined the quantities $\mathrm {i} u^c_{ij} = \mathrm {i} u_{ij} J_{\triangle (\unicode{x2394} )}$ and $\mathrm {i} u^d_{ij} = \mathrm {i} u_{ij} J'_{\triangle (\unicode{x2394} )}$ for $ij \in \triangle (\unicode{x2394} )$. This type of Hamiltonian is solved by taking a Bogoliubov transformation in momentum space. Taking the Fourier transform of (6) and introducing the Nambu spinor as $\Psi^{\rm T}(k) = ( a_1 (k) a_2 (k) \cdots a_6 (k) a^\dagger _1 (-k) \cdots a^\dagger _6 (-k))$, where a_i(k) is the Fourier mode of the fermion creation operator on the ith site in the unit cell, we obtain

$$\begin{equation} {\tilde {\mathcal {H}}} = \sum_{k \in C_+} \Psi^\dagger (k) H(k) \Psi (k), \end{equation} \tag{ 7 }$$

with the identifications

$$\begin{eqnarray} H(k) =\left( \begin{array}{@{}cc@{}} h(k) & \Delta(k)\\ \Delta(k) & -h^{\mathrm{T}}(-k) \end{array}\right),\nonumber \\ &&h_{ij}(k) = \mathrm{i} \left(\tilde{u}_{ij}^c(k) + \tilde{u}_{ij}^d(k)\right) + 2 \hspace{.4mm} J_5 \delta_{ij}, \\ &&\Delta_{ij}(k) = \textstyle\frac{1}{2} \mathrm{i} \left(\tilde{u}^c_{ij}(k) - \tilde{u}^d_{ij}(k)\right),\nonumber \end{eqnarray} \tag{ 8 }$$

where the $\tilde {u}$s are the Fourier modes of the corresponding us defined above, and where C₊ is half of the Brillouin zone (BZ). We only sum over half of the BZ to avoid the double counting introduced by the Nambu spinor. Because H(k) = −H(−k), the eigenvalues of H(k) appear in pairs {E_j(k),−E_j(−k)}, j = 1,2,...,6, and our model is particle–hole symmetric. This is a consequence of our Hamiltonian being written in terms of the Nambu spinor, which satisfies Ψ^†_i(−k) = Ψ_i(k); when we diagonalize H(k), half of the states are redundant.

For general couplings, the ground state may be gapped or gapless. For the gapped phases, the band structure may possess a nontrivial (nonzero) Chern number. We will refer to these nontrivial gapped phases as Chern phases. The physical consequence of a nonzero Chern number is the appearance of chiral gapless edge modes in a system with a boundary. If the system has a finite Chern number, it will exhibit a quantized thermal Hall conductance. We find these signatures of a nontrivial topological phase by calculating the Chern number numerically [29, 43] and diagonalizing a system with boundaries (a 'strip' geometry) to determine the existence of gapless edge modes. This is most conveniently done by taking a Fourier transform to momentum space along the 'length' of the strip geometry [5, 11, 12].

An example band structure is shown in figure 2, which illustrates the dispersion of Majorana excitations $E(\vec k)$ as a function of the two planar momentum components k_x,k_y. The state in figure 2 is gapped about zero energy (below which all states are occupied and above which all states are empty) and has a nontrivial Chern number. There are six bands above and six bands below zero energy, corresponding to the six sites in the unit cell of the ruby lattice. Shown in figure 3 is the spectrum of a strip geometry in the nontrivial phase with Chern number −1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When the ground state is gapless, it generically has a Fermi surface (FS). Shown in figures 4 and 5 are band structures with an FS realized in the ground state. In figure 4 the FS is in the center of the first BZ, while in figure 5 the FS is near the high-symmetry K and K' points of the hexagonal BZ. Note that the FS is 'connected' once the appropriate reciprocal lattice vectors are used to shift them together. At a phase boundary between two gapped phases with different Chern numbers, the ground state is gapless with Dirac nodes at a discrete set of points in the BZ. An example of such a band structure is shown in figure 6, which has four nodes in the first BZ (one at the zone center and three others at the inversion symmetric M points).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In general, the phase diagram is complicated due to the large number of tunable parameters, even when one insists on translational, sixfold rotation and inversion symmetry. We report on a phase diagram resulting from an interesting selection of parameters which reveals the general richness of the model. We explored the phase diagram for the specific case when the couplings $J_{\triangle }=J_{\unicode{x2394} }$ are set equal to a fixed constant J, and the J' couplings are allowed to vary in such a way that their product remains constant; explicitly $J'_{\unicode{x2394} } \times J'_{\triangle }=J^2$. This phase diagram is shown in figure 7 and is exceptionally rich with many gapped and gapless phases.

Zoom In Zoom Out Reset image size

The regions of the phase diagram are classified by whether or not they are gapped or gapless and whether or not the gapless phases possess a nontrivial Chern number (Chern phase). If gapless, the phases are classified according to whether the gap closes at discrete Fermi points in the BZ or along a line or FS. Often the former contains Dirac nodes on the BZ boundary or the zone center where the Fermi point may also be a quadratic band touching point. The Fermi points only occur on critical lines separating two Chern phases or a Chern phase and an FS. We will refer to these critical phases as nodal lines (NL).

The phase diagram in figure 7 contains NLs with up to three Fermi points. Since the NLs are obtained by solving a secular equation of a Bloch Hamiltonian which is analytic, the NLs trace out smooth trajectories in the phase diagram. Moreover, the complicated manner in which they cross yields the many Chern phases in figure 7. Most interestingly, we see a three-node NL which closes in on itself to form an ellipse.

Another interesting region is the high-symmetry line at constant $J'_{\triangle }/J=1$ where $J_{\triangle }=J'_{\triangle }=J_{\unicode{x2394} }=J'_{\unicode{x2394} }=J$. This line contains a gapless FS phase in a line between 0.7 ≲ J₅ ≲ 2.0, and there are four NLs which cross this FS line at its endpoints. Along this line there are no nearby Chern phases in the phase diagram. This FS phase is noncritical with regard to tuning $J'_{\triangle }$ and $J'_{\unicode{x2394} }$, appears to be particularly stable and is the FS phase with the largest volume in the phase diagram.

Finally, we note that in the phase diagram, transitions between Chern phases and FS phases involve critical NLs. This transition will lead to FS pockets forming around the nodes located at the points of high symmetry. Figures 4 and 5 show examples of such pockets. There are also transitions between Chern phases and FS phases which do not involve an NL. In this case, the FS becomes gapped by a pairing (Cooper) instability as the system is driven into the Chern phase.

In section 4, we will study the effects of the inclusion of new interaction terms that spoil the exact solvability of the model. These new interactions also motivate the study of the Ising order in one of the gamma matrices, namely Γ⁵_i.

In this section, we consider the following extension to $\mathcal {H}$:

$$\begin{eqnarray} &&\begin{array}{l} \mathcal{H} = \mathcal{H}_0 + \mathcal{H}_5 + \mathcal{H}_{\mathrm{Ising}}, \\ \mathcal{H}_{\mathrm{Ising}} = - \lambda \displaystyle\sum_{\langle ij \rangle} \Gamma^5_i \Gamma^5_j. \end{array} \end{eqnarray} \tag{ 9 }$$

We limit ourselves to the case where λ > 0. For our analysis it is convenient to take the spin–orbit interpretation of the model where we have four local degrees of freedom with basis states $\{ \left | \uparrow \uparrow \right \rangle , \left | \uparrow \downarrow \right \rangle , \left | \downarrow \uparrow \right \rangle , \left | \downarrow \downarrow \right \rangle \}$. On this basis, Γ⁵_i is trivially diagonal with eigenvalue +1 if the 1/2-spins are aligned and eigenvalue −1 if they are anti-aligned. Thus, we can regard Γ⁵ as an Ising spin, albeit one that is doubly degenerate for each Ising polarization. Then, $\mathcal {H}_{\mathrm {Ising}}$, as its name suggests, describes ferromagnetic nearest-neighbor Ising–Ising interactions in the case where λ > 0. It is also trivially true that $[\mathcal {H}_{\mathrm {Ising}} , \Gamma _i^5 ]=0$ for all sites i, which implies that both operators are simultaneously diagonal in the local spin–orbit basis. Thus in the $\mathcal {H}_0 = 0$ limit, the model is trivially exactly solvable with eigenstates given by tensor products of local Γ⁵_i eigenstates. A typical eigenstate would have the form

$$\begin{eqnarray} & &\left|\psi \right\rangle = \begin{array}{c}\\ \bigotimes\\ \scriptstyle{i\in I}\end{array}\left( \alpha_i \left|\uparrow\uparrow\right\rangle +\beta_i \left|\downarrow\downarrow\right\rangle \right)\begin{array}{c}\\ \bigotimes\\ \scriptstyle{j\in J}\end{array}\left( \gamma_j \left|\uparrow\downarrow\right\rangle +\delta_j \left|\downarrow\uparrow\right\rangle \right),\nonumber\\\\ & &|\alpha_i|^2 + |\beta_i|^2 = |\gamma_j|^2 + |\delta_j|^2 = 1,\nonumber \end{eqnarray} \tag{ 10 }$$

where the N sites of the lattice are partitioned into disjoint sets I and J with N_I and N_J sites each, respectively. The energy of such a state is then given by

$$\begin{eqnarray} &&( \mathcal{H}_5 + \mathcal{H}_{\mathrm{Ising}} )\left|\psi \right\rangle = \left\{ (N_I - N_J )J_5 -\lambda (N_{\mathrm{link}} - L_{\mathrm{dom.}})\right\}\left|\psi \right\rangle, \end{eqnarray}\noindent \tag{ 11 }$$

where N_link is the total number of nearest-neighbor links and L_dom. is the total length of all the boundaries between different 'Ising domains'. The arbitrariness in the parameters α,β,γ,δ on every site leads to a macroscopic degeneracy of the state $\left |\psi \right \rangle$; this stems from the local degeneracy of Γ⁵_i. Also, an exact eigenstate of Γ⁵_i will exhibit no fluctuations in the Γ⁵_i observables. Hence, the ground states of $\mathcal {H}_{\mathrm {Ising}} + \mathcal {H}_5$ are trivial Ising ferromagnets with order in 〈Γ⁵_i〉 = ± 1. The sign of the order parameter will depend on the sign of J₅ if it is nonzero. Otherwise, if J₅ = 0, the ground state will spontaneously break the Ising symmetry in Γ⁵_i.

However, when $\mathcal {H}_0 \neq 0$ the situation is no longer trivial and no longer exactly solvable. In particular, $[\mathcal {H}_0 , \Gamma ^5_i ] \neq 0$ suggests that it would be difficult to diagonalize the Hamiltonian in the spin–orbit basis or any basis where Γ⁵_i is diagonal. Nevertheless, one still has $[W_p,\mathcal {H}_{\mathrm {Ising}}]=0$ for all plaquettes p. This implies that the flux invariants are still good quantum numbers for both the GMM and the $\mathcal {H}_{\mathrm {Ising}}$, and it is sensible to solve for the eigenstates within a given flux sector {W_p}. Writing $\mathcal {H}_{\mathrm {Ising}}$ after performing Majorana transformation makes this explicit,

$$\begin{equation} \tilde{\mathcal{H}}_{\mathrm{Ising}} = -4\lambda \sum_{\langle ij \rangle} a_i^\dagger a_i a_j^\dagger a_j + 2\lambda\sum_i a_i^\dagger a_i - 2N\lambda. \end{equation} \tag{ 12 }$$

Hence, $\tilde {\mathcal {H}}_{\mathrm {Ising}}$ does not couple to the $\mathbb {Z}_2$ gauge field degrees of freedom and the fluxes are conserved. If we now consider again the limit where $\mathcal {H}_0=0$, then $\tilde {\mathcal {H}}$ describes a trivial insulating Hamiltonian with no kinetic terms but with local a-number conservation. One can then write down the exact eigenstates in this representation which is just any occupation state of the complex fermion a_i for all sites i where the flux configuration can remain arbitrary. This gives another explanation for the macroscopic degeneracy of |ψ〉, which in the gauge field picture is due to the many degenerate {W_p} sectors. Equivalently, the flux invariant operators are insensitive to whether or not the local 1/2-spins are aligned or anti-aligned; rather, they operate on degrees of freedom that are orthogonal. We can now give also a physical interpretation of the complex fermion a. Namely, its occupation parity P_i = (2a^†_ia_i − 1) is the Γ⁵_i observable whose eigenvalue determines whether or not the local 1/2-spins are aligned or anti-aligned. Also, the a occupation numbers do not fluctuate in an eigenstate, which is consistent with the nonfluctuating behavior of the Ising order parameter.

Now we return to considering general $\mathcal {H}$. Note that $\mathcal {H}_0$ commutes neither with Γ⁵_i nor with its sum total $\sum _i\Gamma _i^5$. This leads to non-conservation of the total a-particle number in the ground state, and we can associate this nonconservation with the pairing terms in the a-representation of the problem. Nevertheless, the product of parities as defined by $\mathcal {P}=\prod _i P_i = \prod _i \Gamma _i^5$ is conserved since $\mathcal {H}_0,\mathcal {H}_5$ and $\mathcal {H}_{\mathrm {Ising}}$ all individually commute with it. In fact $\mathcal {P}$ implements a global $\mathbb {Z}_2$ gauge transform, a_i← − a_i, a^†_i← − a^†_i on every site i. This is merely another manifestation of the conservation of the total a-number modulo 2 by $\mathcal {H}$.

Next, we will consider the opposite limit where $\mathcal {H}_{\mathrm {Ising}} + \mathcal {H}_5 = 0$ and argue that the exact ground states of $\mathcal {H}_0$ will not possess any Ising order in Γ⁵ by themselves, but that the exact ground states of the more general GMM with J₅ ≠ 0 will. Moreover, in the limit of small J₅, these ground states are adiabatically connected. Given a fixed flux sector {W_p} and given the existence of a nondegenerate ground state of $\tilde {\mathcal {H}}_0$ which we denote by |Ω〉, we will show that 〈Ω|Γ⁵_i|Ω〉 = 0 for all i.

First consider equation (6) in the limit J₅ = 0 where the gauge fields u_ij have been fixed. Define a linear unitary particle–hole conjugation operator $\mathcal {C}$ by its conjugation on a_i and a^†_i and its action on the a-number vacuum |0〉 as follows:

$$\begin{eqnarray} &&\begin{array}{l} \mathcal{C}a_i\mathcal{C}^\dagger = a_i^\dagger, \quad \mathcal{C}a_i^\dagger\mathcal{C}^\dagger = a_i, \\ \mathcal{C}|0\rangle = \mathrm{e}^{\mathrm{i}\theta}a_i^\dagger \ldots a_N^\dagger |0\rangle. \end{array} \end{eqnarray} \tag{ 13 }$$

θ is a phase which we can fix by demanding that $\mathcal {C}^2|0\rangle =|0\rangle$. For N even, this leads to $\theta =\frac {N\pi }{4}$. These relations then totally specify $\mathcal {C}$ and also imply the relations $\mathcal {C}^\dagger =\mathcal {C}$ and $\mathcal {C}^2=1$. Note that $[\tilde {\mathcal {H}}_0,\mathcal {C}] =0$, implying that if |Ω〉 is nondegenerate, then $\mathcal {C}|\Omega \rangle$ can only differ from |Ω〉 by a phase. Thus, $\langle \Omega |a_i^\dagger a_i | \Omega \rangle = \langle \Omega |\mathcal {C}^\dagger a_i a_i^\dagger \mathcal {C}|\Omega \rangle = \langle \Omega |a_ia_i^\dagger |\Omega \rangle$. Using {a_i,a^†_i} = 1, we then conclude that 〈Ω|a^†_ia_i|Ω〉 = 1/2 and 〈Ω|Γ⁵_i|Ω〉 = 0. Since a^†_ia_i is $\mathbb {Z}_2$ gauge invariant, we expect this to hold true even after projection. One could also heuristically argue that because J₅ couples to a^†_ia_i like a chemical potential, the ground state will have a particle–hole symmetry and is the half-filled Fermi sea. However, the presence of pairing terms a_ia_j + a^†_ia^†_j complicates this argument.

For the model that was solved in section 2 where J₅ ≠ 0 generally, 〈Γ⁵_i〉 may take nonzero values in (0,1] when in the ground state. In this more general situation, the presence of $\mathcal {H}_5$ breaks the symmetry under $\mathcal {C}$ conjugation. In addition, since the ground states of the GMM were shown to not undergo a phase transition as J₅ is tuned from 0, we conclude that the exactly solvable ground states of the GMM may exhibit Ising order in the Γ⁵_i spins. That is, the exactly solvable phase and Ising order are not mutually exclusive. Hence the ground state is not a spin liquid in these observables. But fluctuations in Γ⁵_i remain nontrivial in the exactly solvable phase. This poses the interesting question of to what extent an exactly solvable GMM ground state is adiabatically connected to the trivial Ising states such as that described by |ψ〉 above.

In the limit where $\mathcal {H}_0$ is treated as a perturbation to $\mathcal {H}_5+\mathcal {H}_{\mathrm {Ising}}$, the weakly fluctuating Ising order Γ⁵_i will strongly renormalize the J₅ coupling as seen by the a fermions, which can be seen from equation (12). Thus, we can expect this limit to be equivalent to a large J₅ limit and to be described by a low-energy effective Hamiltonian derived from perturbation theory. Schematically, the Hamiltonian can be written in terms of W_p operators [5, 11, 44],

$$\begin{equation} \mathcal{H}_{\mathrm{eff}} = \sum_p \alpha_p W_p + \cdots, \end{equation} \tag{ 14 }$$

where {α_p} are coupling constants and the dots represent higher order terms. Such an effective theory will break the degeneracy among the different flux sectors and favor certain flux sectors over others. However, TRS is still preserved and must be spontaneously broken by the ground state. If J₅ = 0, then the ground state must also spontaneously break the Ising symmetry of Γ⁵.

In the opposite limit where $\mathcal {H}_{\mathrm {Ising}}$ is the perturbation to $\mathcal {H}_0+\mathcal {H}_5$, the situation is more complicated as the ground state of the exactly solvable phase is, in general, nontrivial. We try to address this question in the next section of the paper where $\mathcal {H}_{\mathrm {Ising}}$ is treated as a perturbation to the exactly solvable GMM within the mean-field approximation.

In this section, we study the effects of $\mathcal {H}_{\mathrm {Ising}}$ as a perturbation to the exactly solvable $\mathcal {H}_0+\mathcal {H}_5$ within the mean-field approximation. Performing a mean-field decomposition with respect to the order parameter 〈Γ⁵_i〉 results in the following mean-field Hamiltonian:

$$\begin{equation} \mathcal{H}_{\mathrm{Ising}} \approx \mathcal{H}_{\mathrm{MF}} = - 4 \lambda \sum_{\langle ij \rangle} \langle \Gamma^5_i \rangle (a_i^\dagger a_i -1/2 ) + \lambda \sum_{\langle ij \rangle} \langle \Gamma^5_i \rangle \langle \Gamma^5_j \rangle, \end{equation} \tag{ 15 }$$

where the mean field 〈Γ⁵_i〉 has to be determined self-consistently. Motivated by the fact that the Ising couplings are ferromagnetic, we make the ansatz that 〈Γ⁵_i〉 is uniform across the entire lattice. At the mean-field level, the effect of the interaction $\mathcal {H}_{\mathrm {Ising}}$ can already be seen to renormalize J₅ with J^eff₅ = J₅ + 4λ〈Γ⁵〉.

Within this approximation, we determined the phase diagrams for varying J₅ and λ with fixed J and J' couplings and fixed flux sector. Two phase diagrams are shown in figures 8 and 9. We note that generically the topological phases are stable against weak $\mathcal {H}_{\mathrm {Ising}}$ perturbations. At larger λ, however, the energetics always prefer the trivial gapped phase which has zero Chern number and large 〈Γ⁵〉. Hence, large λ eventually renormalizes J₅ to larger effective values. This is consistent with the large J₅ and large λ trivial phases being adiabatically connected to each other.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Focusing along the line at constant λ = 0 and increasing J₅, we see a series of gapped and gapless phases, which eventually ends with a trivial phase at large J₅. This is as expected from the results of the previous section, but since J₅ acts as an external field, the Ising order is 〈Γ⁵_i〉 ≠ 0 on this line. When λ is tuned into greater values, the volume of these intermediate phases diminishes continuously, eventually giving way to the trivial phase. This confirms that the trivial Ising phase is adiabatically connected to the trivial large J₅ phase.

In figure 8, the gapped phase at J₅ = 0 appears to be the most stable with respect to increasing λ. However, for the couplings shown in figure 9, the J₅ = 0 phase gives way to a trivial gapped phase (Chern number zero) which is not adiabatically connected to the trivial phase at large J₅. This phase undergoes a first-order phase transition to the larger J₅ trivial phase as shown in figure 9. Thus, the Ising interaction may choose to stabilize certain gapped phases over the default J₅ = λ = 0 phase. The stabilized phases are highly dependent on the particular values of the other coupling constants.

At constant J₅ > 0, increasing λ may stabilize certain phases that would typically require larger J₅ values. For example, in figure 8 along the line J₅ = 0.75, the system is driven to an FS phase with increasing λ. Other examples of such transitions driven by the Ising interaction may also been seen in figures 8 and 9. Hence, the Ising perturbation treated at mean field may stabilize certain exotic phases which are otherwise only accessible by intermediate values J₅.

In summary, these results show that with increasing λ, the ground state is eventually adiabatically connected to the trivial Ising ordered state of section 4. However, in general, various phase transitions involving other exotic phases must occur before this happens. Moreover, this depends crucially on the initial (λ = 0) exactly solvable phase.

In this work, we studied an exactly solvable Γ-matrix generalization of Kitaev's original spin-1/2 model on the ruby lattice shown in figure 1. Our model can be interpreted as describing either a spin-3/2 system or a double-orbital spin-1/2 model. We have shown that the ruby lattice, with its many possible tunable parameters, exhibits a rich phase diagram with many exotic phases realized by the same Hamiltonian. We have found gapless phases with FSs and Fermi points, as well as gapped topologically nontrivial phases. Next, we studied the effects of including Ising-like interaction terms that destroy the exact solvability of the model. We analyzed the ordered phase that these terms lead to and argued that the ground state in the large J₅ limit of the exactly solvable model is adiabatically connected to it. We then analyzed in detail the general model near the limit of an exactly solvable phase and treated the Ising terms as perturbations within the mean-field approximation. We derived phase diagrams that show that the Ising interactions will eventually favor the trivial Ising ordered phase of the model over the other more exotic phases. We also confirm that the large J₅ exactly solvable phase is adiabatically connected to the trivial Ising ordered phase.

Although this study has focused on a particular exactly solvable spin model and special perturbations on it, the connections previously drawn to other topological phases indicate that our results are rather broadly applicable [45]. Moreover, many directions for further study in the future suggest themselves. These include the following. (i) Considering anti-ferromagnetic Ising interactions (λ < 0) in order to determine whether the interplay between the exactly solvable GMM and the Ising interactions might be different. (ii) The effect of disorder in interacting many-body quantum systems remains a key open problem. The model we study here with its unusually rich phase diagram provides an opportunity to study disorder effects within a well-controlled model [32]. (iii) Doping magnetic systems is one route to high-temperature superconductivity. To date, no work has been reported on the doping of gamma matrix generalizations of the Kitaev model, and only a few published studies exist on the original model [46]. The richness of the model we study here could shed light on doping magnetic systems in much a more general context and could ultimately help guide the discovery of strongly correlated materials exhibiting high-temperature or unconventional superconductivity.

We gratefully acknowledge funding from ARO grant W911NF-09-1-0527 and NSF grant DMR-0955778.