Exact Chiral Spin Liquids and Mean-Field Perturbations of Gamma Matrix Models on the Ruby Lattice

We theoretically study an exactly solvable Gamma matrix generalization of the Kitaev spin model on the ruby lattice, which is a honeycomb lattice with"expanded"vertices and links. We find this model displays an exceptionally rich phase diagram that includes: (i) gapless phases with stable spin fermi surfaces, (ii) gapless phases with low-energy Dirac cones and quadratic band touching points, and (iii) gapped phases with finite Chern numbers possessing the values {\pm}4,{\pm}3,{\pm}2 and {\pm}1. The model is then generalized to include Ising-like interactions that break the exact solvability of the model in a controlled manner. When these terms are dominant, they lead to a trivial Ising ordered phase which is shown to be adiabatically connected to a large coupling limit of the exactly solvable phase. In the limit when these interactions are weak, we treat them within mean-field theory and present the resulting phase diagrams. We discuss the nature of the transitions between various phases. Our results highlight the richness of possible ground states in closely related magnetic systems.


I. INTRODUCTION
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 in which Kitaev's original version can be generalized, spawning many variants, [8][9][10][11][12][13][14][15][16][17][18][19] and its possession of many non-trivial properties [20][21][22][23][24][25][26][27][28][29] even in the presence of disorder. [30][31][32] Moreover, the exact solution of Kitaev models requires no more effort than solving the problem of non-interacting particles moving in a static background magnetic field. Nonetheless, the exact eigenstates of all Kitaev models are non-trivial entangled manybody 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][34][35][36][37][38][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 non-trivial, non-perturbative 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 on the ruby lattice have yielded a topological insulator, 40 fractional quantum anomalous Hall states, 41 and also topological anyons. 42 Here 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.
The paper is organized as follows. In section II, we introduce the Hamiltonian and discuss its solution by mapping onto a system of non-interacting Majoranas moving in an emergent Z 2 gauge field. In section III, we discuss the various ground state phases that the model realizes and present a phase diagram. In section IV we examine the Ising order in one of the gamma matrix operators which may be induced by including Ising-like interactions. In section V we treat these extra interactions within the mean field approximation and discuss their effects on the phase diagram. We then end with conclusions in section VI.

II. GAMMA-MATRIX MODEL AND HAMILTONIAN
We consider spin-3/2 moments located on the sites of a ruby lattice, shown in Fig.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, . . . 5 11,12,16 . These matrices also have a representation in terms of bilinears of terms of the Γ a i given in either Eq.(1) or Eq.(2). The couplings J , J , J and J are taken to be positive and are link dependent as shown in Fig.1(b). For generic couplings, the Hamiltonian has translational and six-fold 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 timereversal symmetry (TRS), although TRS will be spontaneously broken in the ground state as we will 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 (Ŵ ), four (Ŵ ), and six point (Ŵ ) operators defined on the triangles, squares, and hexagons of the lattice (as their subscripts suggest). Explicitly, they are given byŴ = Γ 12 where the sites are labeled counter-clockwise and for W 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, where the Majorana fermions have the important prop- The Hamiltonian written in terms of these operators then reduces to that of two species of Majorana fermions, c and d, which are non-interacting and move in a background Z 2 gauge field u ij , Here the fields u ij are defined by The u ij are in general quantum fields with eigenvalues ±1 since u 2 ij = 1, hence their identification with a Z 2 gauge theory. Moreover, we can simultaneously diagonalizeH and {u ij } since they commute. We emphasize that the interacting Hamiltonian (3) with the representation (4) is reduced to the effectively non-interacting Hamiltonian (5) because the u ij behave as constants in each flux sector of theŴ ,Ŵ , andŴ .
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 , where the product is over all sites in the lattice. The original Hamiltonian is obtained through H = PHP .
The Z 2 gauge fields define fluxes φ p via exp(iφ p ) ≡ jk∈p iu jk , where jk is taken counterclockwise on each elementary plaquette p. These fluxes are then related to the conservedŴ 's by W p ∝ ij∈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 = ±1's are uniquely specified up to Z 2 gauge transformations, and the Hamil-tonianH is then diagonalized and projected to yield to eigenstates of H. The ground state is then determined by a many-body Majorana wavefunction minimizing the total energy. Due to translational symmetry,H describes a bandstructure 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 time reversal symmetry, 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 time reversal symmetry. 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 which 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 = J = J = J = J and J 5 = 0. The ground state flux with the least energy is depicted by Fig.1(b). Generally the ground state flux will be a function of J 5 , but for simplicity we will not consider additional flux configurations in this work because an additional term (depending onŴ ,Ŵ ,Ŵ ) 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 zero beyond nearest neighbor sites. Hence, the ground state is also a quantum spin-liquid in these observables.
With the gauge sector of H specified, we can solve our system by mapping the Majorana fermions to complex fermions a i , defined by

This puts our Hamiltonian into the form
where we have defined the quantities iu c ij = iu ij J ( ) and iu d ij = iu ij J ( ) for ij ∈ ( ). This type of Hamiltonian is solved by taking a Bogoliubov transfor-mation in momentum space. Taking the Fourier transform of (6) and introducing the Nambu spinor as Ψ † (k) = a 1 (k) a 2 (k) · · · a 6 (k) a † 1 (−k) · · · a † 6 (−k) where a i (k) is the Fourier mode of the fermion creation operator on the ith site in the unit cell, we obtaiñ with the identifications where theũ's are the Fourier modes of the corresponding u's defined above, and where C + is half of the Brillouin zone. We only sum over half of the Brillouin zone to avoid the double counting introduced by the Nambu spinor. Because .., 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.

III. BAND STRUCTURE AND PHASE DIAGRAM
For general couplings the ground state may be gapped or gapless. For the gapped phases, the band structure may possess a non-trivial (non-zero) Chern number. We will refer to these non-trivial gapped phases as Chern phases. The physical consequence of a non-zero 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 Fig.2, which illustrates the dispersion of Majorana excitations E( k) as a function of the two planar momentum components k x , k y . The state in Fig.2 is gapped about zero energy (below which all states are occupied and above which all states are empty) and has a non-trivial Chern number. There are six bands above and six bands below zero energy, corresponding to the 6 sites in the unit cell of the .. There are 12 bands which exhibit a redundant particle-hole symmetry and the energy spectrum is gapped with Chern number 2.

FIG. 3.
Spectrum of the strip geometry of the system as a function of the momentum ky, parallel to the edge. The system has Chern number -1 and couplings {J , J , J , J , J5}={1.0, 0.6, 1.0, 1.67, 0.2}.
Two chiral edge modes that traverse the band gap at E = 0 are clearly visible and each corresponds to a mode on a single edge. ruby lattice. Shown in Fig.3 is the spectrum of a strip geometry in the non-trivial phase with Chern number -1.
When the ground state is gapless, it generically has a Fermi surface. Shown in Figs. 4 and 5 are band structures with a Fermi surface realized in the ground state. In Fig.4 the Fermi surface is in the center of the first Brillouin zone (BZ), while in Fig.5 the Fermi surface is near the high symmetry K and K points of the hexagonal BZ. Note that the Fermi surface 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 Fig.6, which has four nodes in the first BZ (one at the zone center and three others at the inversion symmetric M points).
In general, the phase diagram is complicated due to the large number of tunable parameters, even when one insists on translational, six-fold 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 where the couplings J = J are held equal to a fixed constant J, and the J couplings are allowed to vary in such a way that their product remains a constant; explicitly J × J = J 2 . This phase diagram is shown in Fig.7 and is exceptionally rich with many gapped and gapless phases.
The regions of the phase diagram are classified into whether or not they are gapped or gapless, and whether or not the gapless phases possess a non-trivial 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 Fermi surface (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 a Fermi surface. We will refer to these critical phases as nodal lines (NL).
The phase diagram in Fig.7 contains NLs with up to 3 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 Fig.7. Most interestingly, we see a 3 node NL which closes in on itself to form an ellipse.
Another interesting region is the high symmetry line at constant J /J = 1 where J = J = J = J = J. This line contains an gapless FS phase in a line between 0.7 J 5 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 non-critical in regards to tuning J and J , 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 Fermi surface pockets forming around the nodes located at the points of high symmetry. Figures 4 and 5 are examples of such pockets. There are also transitions between Chern phases and FS phases which do not involve a NL. In this case the Fermi surface becomes gapped by a pairing (Cooper) instability as the system is driven into the Chern phase.
In the next section of this paper, we will study the effects of the inclusion of new interaction terms which spoil the exact solvability of the model. These new interactions will also motivate the study of the Ising order in one the gamma matrices, namely Γ 5 i .

IV. ISING ORDER OF Γ 5
In this section we consider the following extension to H, We also limit ourselves to case where λ > 0. For our analysis it is convenient to take the spin-orbit interpretation of the model where we have 4 local degrees of freedom with basis states {|↑↑ , |↑↓ , |↓↑ , |↓↓ }. In this basis, Γ 5 i is trivially diagonal with eigenvalue +1 if the 1/2−spins are aligned and eigenvalue −1 if they are antialigned. Thus, we can regard Γ 5 as an Ising spin, albeit one that is doubly degenerate for each Ising polarization. Then H Ising , as its name suggests, describes ferromagnetic nearest neighbor Ising-Ising interactions in the case where λ > 0. It is also trivially true that [H Ising , Γ 5 i ] = 0 for all sites i, which implies that both operators are simultaneously diagonal in the local spin-orbit basis. Thus in the H 0 = 0 limit, the model is trivially exactly solvable with eigenstates given by tensor products of local Γ 5 i eigenstates. A typical eigenstate would have the form, 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 where N link is the total number of nearest neighbor links and L dom. is the total length of all the boundaries between the different "Ising domains". The arbitrariness in the parameters α, β, γ, δ on every site leads to a macroscopic degeneracy of the state |ψ ; this stems from the local degeneracy of Γ 5 i . Also, an exact eigenstate of Γ 5 i will exhibit no fluctuations in the Γ 5 i observables. Hence, the ground states of H Ising + H 5 are trivial Ising ferromagnets with order in Γ 5 i = ±1. The sign of the order parameter will depend on the sign of J 5 if it is non-zero. Otherwise, if J 5 = 0, the ground state will spontaneously break the Ising symmetry in Γ Hence,H Ising does not couple to the Z 2 gauge field degrees of freedom and the fluxes are conserved. If we now consider again the limit where H 0 = 0, thenH 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 of 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 also give a physical interpretation to the complex fermion a. Namely, its occupation parity P i = (2a † i a i − 1) is the Γ 5 i observable whose eigenvalue determines whether or not the local 1/2-spin are aligned or anti-aligned. Also, the a occupation numbers do not fluctuate in an eigenstate, which is consistent with the non-fluctuating behavior of the Ising order parameter. Now we return to considering general H. Note that H 0 neither commutes with Γ 5 i nor its sum total i Γ 5 i . This leads to the non-conservation of the total a-particle number in the ground state, and we can associate this non-conservation with the pairing terms in the arepresentation of the problem. Nevertheless, the product of parities as defined by P = i P i = i Γ 5 i is conserved since H 0 , H 5 and H Ising all individually commute with it. In fact P implements a global 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 anumber modulo 2 by H.
Next we will consider the opposite limit where H Ising + H 5 = 0 and argue that the exact ground states of H 0 will not possess any Ising order in Γ 5 by themselves, but that the exact ground states of the more general GMM with J 5 = 0 will. Moreover, in the limit of small J 5 , these ground states are adiabatically connected. Given a fixed flux sector {W p }, and given the existence of a nondegenerate ground state ofH 0 which we denote by |Ω , we will show that Ω|Γ 5 i |Ω = 0 for all i. First consider equation (6) in the limit J 5 = 0 where the gauge fields u ij have been fixed. Define a linear unitary particle-hole conjugation operator C by its by conjugation on a i and a † i and its action on the a-number vacuum |0 as follows: θ is a phase which we can fix by demanding that C 2 |0 = |0 . For N even this leads to θ = N π 4 . These relations then totally specify C and also imply the relations C † = C and C 2 = 1. Note that [H 0 , C] = 0, implying that if |Ω is non-degenerate, then C|Ω can only differ from |Ω by a phase. Thus, Ω|a † i a i |Ω = Ω|C † a i a † i C|Ω = Ω|a i a † i |Ω . Using {a i , a † i } = 1, we then conclude that Ω|a † i a i |Ω = 1/2 and Ω|Γ 5 i |Ω = 0. Since a † i a i is Z 2 gauge invariant, we expect this to hold true even after projection. One could also heuristically argue that because J 5 couples to a † i a 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 i a j + a † i a † j complicates this argument.
For the model that was solved in Sec. II where J 5 = 0 generally, Γ 5 i may take non-zero values in (0, 1] when in the ground state. In this more general situation, the presence of H 5 breaks the symmetry under C conjugation. In addition, since the ground states of the GMM were shown not to undergo a phase transition as J 5 is tuned from zero, we conclude that the exactly solvable ground states of the GMM may exhibit Ising order in the Γ 5 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 Γ 5 i remain non-trivial in the exactly solvable phase. This poses the interesting question as to what extent an exactly solvable GMM ground state is adiabatically connected to the trivial Ising states such as the one described by |ψ above. In the limit where H 0 is treated as a perturbation to H 5 + H Ising , the weakly fluctuating Ising order Γ 5 i will strongly renormalize the J 5 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 5 limit and to be described by a low energy effective Hamiltonian derived from perturbation theory. Schematically, the Hamiltonian can written in terms of W p operators, 5,11,44 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, time reversal symmetry is still preserved and must be spontaneously broken by the ground state. If J 5 = 0, then the ground must also spontaneously break the Ising symmetry of Γ 5 .
In the opposite limit where H Ising is the perturbation to H 0 + H 5 , the situation is more complicated as the ground state of the exactly solvable phase is in general non-trivial. We try to address this question in the next section of the paper where H Ising is treated as a perturbation to the exactly solvable GMM within the mean-field approximation.

V. MEAN FIELD APPROXIMATION OF THE ISING PERTURBATION
In this last section we study the effects of H Ising as a perturbation to the exactly solvable H 0 + H 5 within the mean-field approximation. Performing a mean-field decomposition with respect to the order parameter Γ 5 i results in the following mean-field Hamiltonian, where the mean field Γ 5 i has to be determined selfconsistently. Motivated by the fact that the Ising couplings are ferromagnetic, we make the ansatz that Γ 5 i is uniform across the entire lattice. At the mean-field level, the effect of the interaction H Ising can already be seen to renormalize J 5 with J eff 5 = J 5 + 2λ Γ 5 . Within this approximation, we determined the phase diagrams for varying J 5 and λ with fixed J and J couplings and fixed flux sector. Two phase diagrams are shown in Fig.8 and Fig.9. We note that generically the topological phases are stable against weak H Ising perturbations. At larger λ, however, the energetics always prefer the trivial gapped phase which has zero Chern number and large Γ 5 . Hence, large λ eventually renormalizes J 5 to larger effective values. This is consistent with the large J 5 and large λ trivial phases being adiabatically connected to each other.
Focusing along the line at constant λ = 0 and increasing J 5 , we see a series of gapped and gapless phases which eventually ends with a trivial phase at large J 5 . This is expected from the results of the previous section, but since J 5 acts like an external field, the Ising order is Γ 5 i = 0 on this line. When λ is tuned to 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 5 phase.
In Fig.8, the gapped phase at J 5 = 0 appears to be the most stable with respect to increasing λ. However, for the couplings shown in Fig.9, the J 5 = 0 phase gives way to a trivial gapped phase (Chern number zero) which is not adiabatically connected to the trivial phase at large J 5 . This phase undergoes a first order phase transition to the larger J 5 trivial phase as shown in Fig.9. Thus, the Ising interaction may choose to stabilize certain gapped phases over the default J 5 = λ = 0 phase. The stabilized phases are highly dependent on the particular values of the other coupling constants.
At constant J 5 > 0, increasing λ may stabilize certain phases which would typically require larger J 5 values. For example, in Fig.8 along the line J 5 = 0.75, the system is driven to a Fermi-surface phase with increasing λ.
Other examples of such transitions driven by the Ising interaction may also been seen in Figs. 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 5 .
In summary, these results show that with increasing λ, the ground state is eventually adiabatically connected to the trivial Ising ordered state of Sec. IV. 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.

VI. CONCLUSIONS
In this work we studied an exactly solvable Γ-matrix generalization of Kitaev's original spin-1/2 model on the ruby lattice shown in Fig.1. Our model can be interpreted as describing either a spin-3/2 system or a doubleorbital 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 same Hamiltonian. We have found gapless phases with Fermi surfaces and Fermi points, as well as gapped topologically non-trivial 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 5 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 treat the Ising terms as perturbations within the mean-field approximation. We derived phase diagrams which 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 5 exactly solvable phase is adiabatically connected to the trivial Ising ordered phase.
While this study has focused on a particular exactly solvable spin model and special perturbations on it, connections previously drawn to other topological phases indicate that its results are rather broadly applicable. 45 Moreover, a number of future directions for further study suggest themselves: (i) Considering anti-ferromagnetic Ising interactions (λ < 0) to determine if the interplay between the exactly solvable GMM and the Ising in-teractions 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 studies exist on the original model. 46 The richness of the model we study here could shed light on doping mag-netic systems in a much more general context and could ultimately help guide the discovery of strongly correlated materials exhibiting high temperature or unconventional superconductivity.