Vison Crystals in an Extended Kitaev Model on the Honeycomb Lattice
Shang-Shun Zhang, Zhentao Wang, Gábor B. Halász, Cristian D. Batista
VVison Crystals in an Extended Kitaev Model on the Honeycomb Lattice
Shang-Shun Zhang, Zhentao Wang, G´abor B. Hal´asz, and Cristian D. Batista
1, 3 Department of Physics and Astronomy, The University of Tennessee, Knoxville, Tennessee 37996, USA Materials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA Neutron Scattering Division and Shull-Wollan Center,Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA
We introduce an extension of the Kitaev honeycomb model by including four-spin interactions thatpreserve the local gauge structure and hence the integrability of the original model. The extendedmodel has a rich phase diagram containing five distinct vison crystals, as well as a symmetric π -fluxspin liquid with a Fermi surface of Majorana fermions and a sequence of Lifshitz transitions. Wediscuss possible experimental signatures and, in particular, present finite-temperature Monte Carlocalculations of the specific heat and the static vison structure factor. We argue that our extendedmodel emerges naturally from generic perturbations to the Kitaev honeycomb model. Introduction.
The famous Kitaev model on the honey-comb lattice [1] is an exactly solvable yet experimentallyrealistic model of a quantum spin liquid. In contrast tomore conventional magnetic phases, quantum spin liquidsretain extensive (quantum) fluctuations all the way downto zero temperature [2], where the spins appear to frac-tionalize into deconfined “spinon” quasiparticles coupledto appropriate gauge fields [3].The Kitaev model is approximately realized in a fam-ily of strongly spin-orbit-coupled honeycomb materials,where its anisotropic spin interactions emerge betweeneffective J = 1 / t g orbitalsof 4 d or 5 d ions [4–7]. To determine the most accuratemicroscopic spin models for these Kitaev materials, in-cluding (Na,Li) IrO [8–16] and α -RuCl [17–30], variousextensions of the Kitaev model have been considered andanalyzed with a wide range of techniques [31–51]. Whilethese models are experimentally realistic and have richphase diagrams in the classical limit, it is challenging toidentify and characterize quantum phases in them. Fora start, the honeycomb lattice may harbor many differ-ent quantum spin liquids [52, 53], and the Kitaev spinliquid, captured by the Kitaev model, is only one amongthese many candidates. In addition, a quantum spin liq-uid may also remain “hidden” by appearing on top ofclassical symmetry-breaking order [54].From a more phenomenological point of view, the low-energy physics of the Kitaev spin liquid is described byMajorana fermions (spinons) with Dirac nodes, coupledto an emergent Z gauge field [1]. At each plaquette ofthe honeycomb lattice, the Z gauge field may form a π flux, corresponding to a “vison” excitation. In turn,the presence of such a vison affects the kinetic energyof the spinons via the Berry phase π picked up by eachspinon moving around it. For the pure Kitaev model,the spinons are governed by a nearest-neighbor hoppingproblem (cf. electrons in graphene) and, due to the lack offrustration, the ground state has no visons at any plaque-ttes [1, 55]. However, if the hopping problem is frustratedby competing hopping amplitudes, the presence of a vi-son may reduce the frustration and thus lower the kinetic (b) (c) il i jk l (a) x yz AB j (cid:31) k (cid:31) jkp FIG. 1. Extended Kitaev model. (a) Honeycomb lattice withtwo sublattices A and B (black and white dots), three bondtypes x , y , and z (red, green, and blue bonds), and the site-labeling convention around a plaquette p . (b)-(c) Representa-tive (orange) paths (cid:104) ijkl (cid:105) yzx (b) and (cid:104) ijkl (cid:105) yzy (c) associatedwith the K and K (cid:48) terms in Eq. (2), respectively; four-spininteractions along such paths give rise to Majorana hoppingfrom any site i to all its third neighbors [59], as indicated bythe dashed arrows. For the path (cid:104) ijkl (cid:105) yzx (b), the symmetry-related path (cid:104) ij (cid:48) k (cid:48) l (cid:105) xzy is marked by blue. energy of the spinons. Such a frustration in the hoppingamplitudes is known to stabilize crystals of topologicalsolitons, such as baby skyrmions or merons, in itinerantmagnets [56–58], and one may thus expect it to stabilizeanalogous vison crystals in the Kitaev spin liquid.In this Letter, we extend the Kitaev model by includ-ing four-spin interactions that preserve the exact solutionof the model and emerge naturally from generic pertur-bations. By introducing frustrated further-neighbor hop-ping for the Majorana fermions, these additional inter-actions stabilize a rich variety of vison crystals, as wellas a symmetric π -flux spin liquid with a vison at everyplaquette. Interestingly, the π -flux spin liquid exhibits aFermi surface of Majorana fermions undergoing two sub-sequent Lifshitz transitions. On a technical level, we firstuse a simple variational treatment to compute the zero-temperature phase diagram of our extended model. Thevalidity of this approach is then confirmed by unbiasedMonte Carlo (MC) simulations that also reveal the finitemelting temperatures of the vison crystals. Model.
We consider a generalized Kitaev Hamiltonian a r X i v : . [ c ond - m a t . s t r- e l ] J u l on the honeycomb lattice: H = H K + H K , (1)where H K = − K (cid:80) (cid:104) ij (cid:105) α σ αi σ αj is the usual [1] isotropicKitaev Hamiltonian with ferromagnetic ( K >
0) Isinginteractions between the spin components σ α along each α = { x, y, z } bond (cid:104) ij (cid:105) α [see Fig. 1(a)], and H K = K (cid:88) (cid:104) ijkl (cid:105) αβγ σ αi σ γj σ αk σ γl − K (cid:48) (cid:88) (cid:104) ijkl (cid:105) αβα σ αi σ γj σ γk σ αl , (2)where ( αβγ ) is a permutation of ( xyz ) in each term, and (cid:104) ijkl (cid:105) αβγ is a path of length 3 consisting of bonds (cid:104) ij (cid:105) α , (cid:104) jk (cid:105) β , and (cid:104) kl (cid:105) γ . Each term in H K is the product of thethree terms in H K that correspond to the three bondsalong the appropriate path. Different K and K (cid:48) termsare related by space-group symmetries, simultaneouslytransforming the lattice and the spins; particular exam-ples of their respective paths, with ( αβγ ) = ( yzx ), aredepicted in Figs. 1(b) and 1(c). We remark that, for eachpath (cid:104) ijkl (cid:105) αβγ going around one “half” of a hexagon,connecting opposite vertices i and l , there is a symmetry-related path (cid:104) lk (cid:48) j (cid:48) i (cid:105) αβγ = (cid:104) ij (cid:48) k (cid:48) l (cid:105) γβα going around theother “half” of the hexagon [see Fig. 1(b)].Importantly, the exact solution of H K [1] is preservedby the additional terms in Eq. (2). Indeed, since H com-mutes with the flux operator W p = σ x σ y σ z σ x σ y σ z ateach plaquette p [see Fig. 1(a)], one can identify static Z flux or “vison” degrees of freedom at these plaquettes,each being present (absent) if the corresponding W p takeseigenvalue − σ αj = ib αj c j , the Hamiltonian takes the form [60] H = iK (cid:88) (cid:104) ij (cid:105) α u αij c i c j + iK (cid:88) (cid:104) ijkl (cid:105) αβγ u αij u βkj u γkl c i c l + iK (cid:48) (cid:88) (cid:104) ijkl (cid:105) αβα u αij u βkj u αkl c i c l , (3)where u αij = − u αji ≡ ib αi b αj is a Z gauge field along the α bond (cid:104) ij (cid:105) α . Since these gauge fields are conserved quan-tities, u αij = ±
1, providing a redundant description of theconserved gauge fluxes, W p = u z u x u y u z u x u y = ± c i , thusgiving rise to free fermion (“spinon”) excitations after astraightforward diagonalization [60]. From the perspec-tive of the Majorana fermions, the K terms describefirst-neighbor hopping, while the additional K and K (cid:48) terms describe third-neighbor hopping [59].In analogy with how three-spin interactions may be ob-tained from a Zeeman field [1], the four-spin interactionsin Eq. (2) can in principle be generated by a perturbativetreatment of Heisenberg and/or symmetric off-diagonal(Γ) interactions on top of the pure Kitaev model. Tak-ing a more universal approach and considering Eq. (3) asan effective low-energy theory for the Majorana fermions Zero flux(Dirac) π flux (Fermi) π flux (Fermi) K (cid:31) / K K /K π flux (Dirac) 2/3 flux (gapped)1/3 flux (gapped)2/3 flux (Dirac) FIG. 2. Phase diagram of the extended Kitaev model. Fluxconfigurations of distinct vison crystals (colored phases) aredepicted in separate panels; the presence (absence) of a fluxis marked at each plaquette by gray (white) filling. [61], we know that generic time-reversal-symmetric per-turbations to H K must generate all Majorana terms thatare consistent with the projective symmetries of the Ki-taev spin liquid [53]. Given that all interaction terms areirrelevant and second-neighbor hopping terms are forbid-den by time reversal, Eq. (3) is the most natural effectivetheory beyond the pure Kitaev model. Phase diagram.
The ground state of H K belongs tothe zero-flux sector, characterized by W p = +1 for all p [1, 55]. In the presence of the additional interactions,however, the ground state may belong to a wide rangeof different flux sectors, as shown by the T = 0 phasediagram in Fig. 2. This phase diagram is obtained froma simple variational analysis, by comparing the energiesof the seven flux sectors appearing in the diagram onfinite lattices of 48 ×
48 unit cells [62]. Furthermore, itis fully consistent with unbiased finite-temperature MCsimulations, discussed in a later section [63].We first concentrate on the two fully symmetric non-crystal phases occupying most of the phase diagram: thezero-flux phase, which has no fluxes at any plaquettes,and the π -flux phase, which has a Z flux at each pla-quette. For K = K (cid:48) = 0, the creation of each Z fluxwith W p = − ≈ . K , and theground state thus belongs to the zero-flux sector. For K /K >
0, the K and K terms in Eq. (3) give rise toa frustrated Majorana hopping and hence an increase inthe ground-state energy. However, due to the two pathsbetween any two opposite sites i and l around a plaquette p [see Fig. 1(b)], there are two equivalent hopping terms ∝ iK c i c l in Eq. (3), which interfere constructively for W p = +1 and destructively for W p = −
1. Consequently,as K /K is increased, fluxes are effective in relieving FIG. 3. Majorana nodal structures (dark blue) in the various π -flux phases: the Dirac phase (a), the first Fermi phase (b),the Lifshitz transition between the two Fermi phases (c), andthe second Fermi phase (d). In the presence of generic further-neighbor Majorana hopping terms [59], each Fermi surface isgapped out into six Dirac points (red circles). frustration from the Majorana hopping and thus becomeenergetically favorable. Since the effective interaction be-tween nearby fluxes is attractive for small K (cid:48) /K [1], thecorresponding phase transition between the zero-flux andthe π -flux phases is strongly first order.Increasing K (cid:48) /K , one can modify this interaction andstabilize various intermediate phases with nontrivial fluxconfigurations. Indeed, there are five distinct translation-symmetry-breaking vison-crystal phases in Fig. 2, withtheir ordering wave vectors Q corresponding to either theK point or the M point(s) of the Brillouin zone (BZ). Thetwo Q = Q K crystals have supercells of three plaquettes,containing one vison (“1 / / Q = Q M crystals can exhibit single- Q or multi- Q ordering. The single- Q crystal is a stripy con-figuration, corresponding to a supercell of two plaquettescontaining one vison (“1 / Q crystals have supercells of four plaquettes, con-taining one vison (“1 / / Majorana problems.
For the different ground-state fluxsectors discussed above, distinct configurations of thegauge fields u αij = ± / / / / / π -flux phase. This phase isamenable to a full analytic understanding as, due to theperfect cancelation of all K terms in Eq. (3), the Majo-rana problem has only one dimensionless parameter ratio κ ≡ K (cid:48) /K . With a simple calculation [60], we find thatthere are in fact three distinct π -flux phases characterizedby different Majorana nodal structures.In particular, there is a π -flux phase where the Majo-rana fermions are gapless at Dirac points only, and an-other two π -flux phases where these Dirac points coexistwith Fermi surfaces (i.e., nodal lines) of distinct topolo-gies (see Fig. 3). The dashed lines in Fig. 2 indicate twosubsequent Lifshitz transitions [64] separating these threephases as a function of increasing κ . For κ < /
5, theonly nodal structures are Dirac points. At the first Lif-shitz transition, κ = 1 /
5, small pockets of Fermi surfacesappear around these Dirac points and gradually expandas κ is further increased. At the second Lifshitz transi-tion, κ = ( √ − / ≈ . κ .Such a coexistence of Dirac points and Fermi surfacesis rather surprising and is not expected to be stable. In-stead, due to the nature of the time-reversal and particle-hole symmetries in the Majorana problem [65], one wouldanticipate only Dirac points to be generically present, asin all the other phases of Fig. 2. Indeed, we find thatthe Fermi surfaces exist due to the particular simplicityof the problem up to third-neighbor hopping terms [60]and that each Fermi surface is gapped out into six Diracpoints (see Fig. 3) when generic fifth-neighbor hoppingterms [59], respecting the projective symmetries of thesystem, are included in Eq. (3). However, assuming thatsuch terms are small enough, approximate Fermi surfacesare still expected to be observable in experiments.
Experimental signatures.
The phase diagram in Fig. 2contains a rich variety of phases with all possible Majo-rana nodal structures in two dimensions, including Fermisurfaces, Dirac points, and fully gapped scenarios. Due totheir distinct low-energy physics, these phases are char-acterized by different experimental signatures. First, weexpect the low-temperature specific heat to behave as C ∝ T for Fermi phases, C ∝ T for Dirac phases, and C ∝ e − ∆ v /T for fully gapped phases, where the activated FIG. 4. Dynamical spin structure factor S zz ( q , ω ) [60] for the1/3 flux crystal (a) and the 3/4 flux crystal (b) along the pathM-Γ-K-M in the Brillouin zone [see inset of panel (a)] via thesingle-particle approximation of Ref. 69. behavior should be controlled by the vison gap ∆ v as itis actually smaller than the Majorana gap. Second, thevarious Majorana nodal structures may be distinguishedby their low-energy fingerprints in spectroscopic probes,such as resonant inelastic x-ray scattering [66, 67]. Third,the Majorana Fermi surface in the π -flux phase leads toimpurity-induced Friedel oscillations in the magnetic en-ergy density [60]. In turn, such magnetic Friedel oscilla-tions should be measurable with nuclear magnetic reso-nance (NMR) as they induce an oscillatory bond-lengthmodulation via magnetostriction.For the vison-crystal phases in Fig. 2, the spontaneousbreaking of translation symmetry leads to further exper-imental signatures. First of all, due to magnetostriction,each vison crystal generates a characteristic bond-lengthmodulation throughout the lattice, which can be pickedup with NMR or elastic x-ray scattering. Moreover, theenlarged unit cell results in a larger number of distinctbands for the Majorana fermions and therefore, in con-trast to the pure Kitaev model [68, 69], the dynamicalspin structure factor [60], directly measurable by inelas-tic neutron scattering, has multiple peaks as a functionof energy (see Fig. 4). Finally, unlike the fully symmetricphases, each vison-crystal phase has a finite-temperaturephase transition at a critical temperature T c . Monte Carlo simulations.
To verify the phase diagramin Fig. 2 and to extract the melting temperatures T c ofthe vison crystals, we perform MC simulations of H basedon a Metropolis algorithm to update the “classical” Z fields { u αij = ± } . The energy of each field configura-tion is computed by diagonalizing the quadratic Majo-rana Hamiltonian in Eq. (3) [70, 71] on L × L latticeswith L = { , , } [60]. For each temperature, a sin-gle run contains 10000 MC sweeps for equilibration andanother 20000 MC sweeps for measurement [72].Figure 5 shows our results for the heat capacity C ( T )and the static vison structure factor, ρ v ( k ) = 1 L (cid:88) p,p (cid:48) e i k · ( X p − X p (cid:48) ) (cid:104) W p W p (cid:48) (cid:105) , (4)for representative parameters of four different vison crys-tals, where X p is the position of plaquette p , and k is the . . . . / / − − − (c)2 / − − − (d)3 / − − − − − − − − − − − C / NC / N T / K L = L = L = T / K ρ v ( K ) / N ρ v ( M ) / N ρ v ( K ) / N ρ v ( M ) / N FIG. 5. Temperature dependence of the specific heat and theappropriate static vison structure factor for (a) the 1 / K = 0 . K and K (cid:48) = 0 . K , (b) the 1 / K = 0 . K and K (cid:48) = 0 . K , (c) the2 / K = 0 . K and K (cid:48) = 0 . K , and(d) the 3 / K = 0 . K and K (cid:48) = 0 . K on L × L lattices ( L = 6 , ,
18) containing N = 2 L sites. ordering wave vector of each vison crystal, correspondingto either the K or the M point of the BZ. We first observethat, as for the pure Kitaev model, C ( T ) exhibits both ahigh- and a low-temperature peak, which correspond tospinon and vison excitations, respectively [70, 71]. How-ever, the low-temperature peak signals the onset of vison-crystal ordering at T = T c , as confirmed by the sharpgrowth of the corresponding Bragg peak in ρ v ( k ). Whilethe three lattice sizes L = { , , } do not facilitate arigorous finite-size scaling analysis, the results in Fig. 5suggest a first-order crystallization transition for all visoncrystals, except for the 1 / / C ( T ) /L shouldbe ∝ L α/ν with critical exponents α = 1 / ν = 5 /
6, and α/ν = 2 / T c ∼ − K . Discussion.
By considering a natural extension of thehoneycomb Kitaev model, we have found a rich spec-trum of novel spin-liquid phases that are not adiabati-cally connected to the original Kitaev model, includinga fully symmetric π -flux spin liquid, and five distinctsymmetry-breaking spin liquids with various degrees ofvison crystallization. In the future, it would be interest-ing to study how an external magnetic field affects ourspin liquids. For the Dirac phases, it may generate non-Abelian gapped spin liquids with distinct Chern numbersof the Majorana fermions [1]. For the gapped phases, itmay lead to nontrivial finite-field phase transitions be-tween topologically distinct spin liquids.We thank Arnab Banerjee, Hiroaki Ishizuka, and Jo-hannes Knolle for useful comments on the manuscript.S.-S.Z., Z.W., and C.D.B. are supported by funding fromthe Lincoln Chair of Excellence in Physics. The work ofG.B.H. at ORNL was supported by Laboratory Direc-tor’s Research and Development funds. This researchused resources of the Oak Ridge Leadership ComputingFacility, which is a DOE Office of Science User Facilitysupported under Contract DE-AC05-00OR22725. [1] A. Kitaev, Ann. Phys. , 2 (2006).[2] L. Balents, Nature , 199 (2010).[3] L. Savary and L. Balents, Rep. Prog. Phys. , 016502(2016).[4] G. Jackeli and G. Khaliullin, Phys. Rev. Lett. ,017205 (2009).[5] J. G. Rau, E. K.-H. Lee, and H.-Y. Kee, Annu. Rev.Condens. Matter Phys. , 195 (2016).[6] S. Trebst, ArXiv e-prints (2017), arXiv:1701.07056[cond-mat.str-el].[7] M. Hermanns, I. Kimchi, and J. Knolle, Annu. Rev.Condens. Matter Phys. , 17 (2018).[8] Y. Singh and P. Gegenwart, Phys. Rev. B , 064412(2010).[9] X. Liu, T. Berlijn, W.-G. Yin, W. Ku, A. Tsvelik, Y.-J.Kim, H. Gretarsson, Y. Singh, P. Gegenwart, and J. P.Hill, Phys. Rev. B , 220403 (2011).[10] Y. Singh, S. Manni, J. Reuther, T. Berlijn, R. Thomale,W. Ku, S. Trebst, and P. Gegenwart, Phys. Rev. Lett. , 127203 (2012).[11] S. K. Choi, R. Coldea, A. N. Kolmogorov, T. Lancaster,I. I. Mazin, S. J. Blundell, P. G. Radaelli, Y. Singh,P. Gegenwart, K. R. Choi, S.-W. Cheong, P. J. Baker,C. Stock, and J. Taylor, Phys. Rev. Lett. , 127204(2012).[12] F. Ye, S. Chi, H. Cao, B. C. Chakoumakos, J. A.Fernandez-Baca, R. Custelcean, T. F. Qi, O. B. Korneta,and G. Cao, Phys. Rev. B , 180403 (2012).[13] R. Comin, G. Levy, B. Ludbrook, Z.-H. Zhu, C. N. Veen-stra, J. A. Rosen, Y. Singh, P. Gegenwart, D. Stricker,J. N. Hancock, D. van der Marel, I. S. Elfimov, andA. Damascelli, Phys. Rev. Lett. , 266406 (2012).[14] S. Hwan Chun, J.-W. Kim, J. Kim, H. Zheng, C. C.Stoumpos, C. D. Malliakas, J. F. Mitchell, K. Mehlawat,Y. Singh, Y. Choi, T. Gog, A. Al-Zein, M. M. Sala,M. Krisch, J. Chaloupka, G. Jackeli, G. Khaliullin, andB. J. Kim, Nat. Phys. , 462 (2015).[15] S. C. Williams, R. D. Johnson, F. Freund, S. Choi,A. Jesche, I. Kimchi, S. Manni, A. Bombardi, P. Manuel,P. Gegenwart, and R. Coldea, Phys. Rev. B , 195158(2016).[16] K. Kitagawa, T. Takayama, Y. Matsumoto, A. Kato,R. Takano, Y. Kishimoto, R. Dinnebier, G. Jackeli, andH. Takagi, Nature , 341 (2018).[17] K. W. Plumb, J. P. Clancy, L. J. Sandilands, V. V.Shankar, Y. F. Hu, K. S. Burch, H.-Y. Kee, and Y.-J. Kim, Phys. Rev. B , 041112 (2014).[18] L. J. Sandilands, Y. Tian, K. W. Plumb, Y.-J. Kim, andK. S. Burch, Phys. Rev. Lett. , 147201 (2015). [19] J. A. Sears, M. Songvilay, K. W. Plumb, J. P. Clancy,Y. Qiu, Y. Zhao, D. Parshall, and Y.-J. Kim, Phys.Rev. B , 144420 (2015).[20] M. Majumder, M. Schmidt, H. Rosner, A. A. Tsirlin,H. Yasuoka, and M. Baenitz, Phys. Rev. B , 180401(2015).[21] R. D. Johnson, S. C. Williams, A. A. Haghighirad, J. Sin-gleton, V. Zapf, P. Manuel, I. I. Mazin, Y. Li, H. O.Jeschke, R. Valent´ı, and R. Coldea, Phys. Rev. B ,235119 (2015).[22] L. J. Sandilands, Y. Tian, A. A. Reijnders, H.-S. Kim,K. W. Plumb, Y.-J. Kim, H.-Y. Kee, and K. S. Burch,Phys. Rev. B , 075144 (2016).[23] A. Banerjee, C. A. Bridges, J.-Q. Yan, A. A. Aczel, L. Li,M. B. Stone, G. E. Granroth, M. D. Lumsden, Y. Yiu,J. Knolle, S. Bhattacharjee, D. L. Kovrizhin, R. Moess-ner, D. A. Tennant, M. D. G., and S. E. Nagler, Naturematerials (2016), 10.1038/nmat4604.[24] J. A. Sears, Y. Zhao, Z. Xu, J. W. Lynn, and Y.-J. Kim,Phys. Rev. B , 180411 (2017).[25] A. Banerjee, J. Yan, J. Knolle, C. A. Bridges, M. B.Stone, M. D. Lumsden, D. G. Mandrus, D. A. Tennant,R. Moessner, and S. E. Nagler, Science , 1055 (2017).[26] S.-H. Baek, S.-H. Do, K.-Y. Choi, Y. S. Kwon, A. U. B.Wolter, S. Nishimoto, J. van den Brink, and B. B¨uchner,Phys. Rev. Lett. , 037201 (2017).[27] S.-H. Do, S.-Y. Park, J. Yoshitake, J. Nasu, Y. Motome,Y. S. Kwon, D. T. Adroja, D. J. Voneshen, K. Kim, T.-H. Jang, J.-H. Park, K.-Y. Choi, and S. Ji, Nat. Phys. , 1079 (2017).[28] A. Banerjee, P. Lampen-Kelley, J. Knolle, C. Balz, A. A.Aczel, B. Winn, Y. Liu, D. Pajerowski, J. Yan, C. A.Bridges, A. T. Savici, B. C. Chakoumakos, M. D. Lums-den, D. A. Tennant, R. Moessner, D. G. Mandrus,and S. E. Nagler, npj Quantum Materials , 8 (2018),arXiv:1706.07003 [cond-mat.mtrl-sci].[29] R. Hentrich, A. U. B. Wolter, X. Zotos, W. Brenig,D. Nowak, A. Isaeva, T. Doert, A. Banerjee, P. Lampen-Kelley, D. G. Mandrus, S. E. Nagler, J. Sears, Y.-J. Kim,B. B¨uchner, and C. Hess, Phys. Rev. Lett. , 117204(2018).[30] Y. Kasahara, T. Ohnishi, Y. Mizukami, O. Tanaka,S. Ma, K. Sugii, N. Kurita, H. Tanaka, J. Nasu, Y. Mo-tome, T. Shibauchi, and Y. Matsuda, Nature , 227(2018).[31] J. Chaloupka, G. Jackeli, and G. Khaliullin, Phys. Rev.Lett. , 027204 (2010).[32] H.-C. Jiang, Z.-C. Gu, X.-L. Qi, and S. Trebst, Phys.Rev. B , 245104 (2011).[33] J. Reuther, R. Thomale, and S. Trebst, Phys. Rev. B , 100406 (2011).[34] C. C. Price and N. B. Perkins, Phys. Rev. Lett. ,187201 (2012).[35] J. G. Rau, E. K.-H. Lee, and H.-Y. Kee, Phys. Rev. Lett. , 077204 (2014).[36] Y. Yamaji, Y. Nomura, M. Kurita, R. Arita, andM. Imada, Phys. Rev. Lett. , 107201 (2014).[37] Y. Sizyuk, C. Price, P. W¨olfle, and N. B. Perkins, Phys.Rev. B , 155126 (2014).[38] E. Sela, H.-C. Jiang, M. H. Gerlach, and S. Trebst, Phys-ical Review B , 035113 (2014).[39] I. Rousochatzakis, J. Reuther, R. Thomale, S. Rachel,and N. B. Perkins, Phys. Rev. X , 041035 (2015).[40] H.-S. Kim, V. V., A. Catuneanu, and H.-Y. Kee, Phys. Rev. B , 241110 (2015).[41] R. Yadav, N. A. Bogdanov, V. M. Katukuri, S. Nishi-moto, J. Van Den Brink, and L. Hozoi, Scientific reports , 37925 (2016).[42] H.-S. Kim and H.-Y. Kee, Physical Review B , 155143(2016).[43] S. M. Winter, Y. Li, H. O. Jeschke, and R. Valent´ı, Phys.Rev. B , 214431 (2016).[44] B. H. Kim, T. Shirakawa, and S. Yunoki, Physical reviewletters , 187201 (2016).[45] L. Janssen, E. C. Andrade, and M. Vojta, Physical Re-view Letters , 277202 (2016).[46] Y. S. Hou, H. J. Xiang, and X. G. Gong, Phys. Rev. B , 054410 (2017).[47] S. M. Winter, A. A. Tsirlin, M. Daghofer, J. van denBrink, Y. Singh, P. Gegenwart, and R. Valent´ı, Journalof Physics: Condensed Matter , 493002 (2017).[48] A. M. Samarakoon, A. Banerjee, S.-S. Zhang, Y. Kamiya,S. E. Nagler, D. A. Tennant, S.-H. Lee, and C. D.Batista, Phys. Rev. B , 134408 (2017).[49] K. Ran, J. Wang, W. Wang, Z.-Y. Dong, X. Ren, S. Bao,S. Li, Z. Ma, Y. Gan, Y. Zhang, et al. , Physical ReviewLetters , 107203 (2017).[50] A. M. Samarakoon, G. Wachtel, Y. Yamaji, D. A. Ten-nant, C. D. Batista, and Y. B. Kim, Phys. Rev. B ,045121 (2018).[51] J. S. Gordon, A. Catuneanu, E. S. Sørensen, andH.-Y. Kee, arXiv e-prints , arXiv:1901.09943 (2019),arXiv:1901.09943 [cond-mat.str-el].[52] Y.-M. Lu and Y. Ran, Phys. Rev. B , 024420 (2011).[53] Y.-Z. You, I. Kimchi, and A. Vishwanath, Phys. Rev. B , 085145 (2012).[54] L. Savary and L. Balents, Phys. Rev. Lett. , 037202(2012).[55] E. H. Lieb, Phys. Rev. Lett. , 2158 (1994).[56] R. Ozawa, S. Hayami, K. Barros, G.-W. Chern, Y. Mo-tome, and C. D. Batista, J. Phys. Soc. Jpn. , 103703(2016).[57] C. D. Batista, S.-Z. Lin, S. Hayami, and Y. Kamiya,Reports on Progress in Physics , 084504 (2016).[58] R. Ozawa, S. Hayami, and Y. Motome, Phys. Rev. Lett. , 147205 (2017).[59] In this Letter, “ n -th-neighbor” means that the shortestpath connecting the two sites consists of n bonds.[60] See Supplemental Material for extended descriptions ofthe quadratic Majorana problems and the correspondingnodal structures in the various phases, for detailed resultson the magnetic Friedel oscillations and the dynamicalspin structure factor, as well as for implementation de-tails of the Monte Carlo simulations.[61] X.-Y. Song, Y.-Z. You, and L. Balents, Phys. Rev. Lett. , 037209 (2016).[62] We verified that fluctuations in the phase boundaries dueto finite-size effects become negligibly small for latticeslarger than 36 ×
36 unit cells.[63] We verified this statement by running unbiased MC sim-ulations for multiple randomly chosen points within eachphase on finite lattices of 12 ×
12 unit cells.[64] G. E. Volovik,
The Universe in a Helium Droplet (OxfordUniversity Press, New York, USA, 2003).[65] M. Hermanns and S. Trebst, Phys. Rev. B , 235102(2014).[66] G. B. Hal´asz, N. B. Perkins, and J. van den Brink, Phys.Rev. Lett. , 127203 (2016). [67] G. B. Hal´asz, B. Perreault, and N. B. Perkins, Phys.Rev. Lett. , 097202 (2017).[68] J. Knolle, D. L. Kovrizhin, J. T. Chalker, and R. Moess-ner, Phys. Rev. Lett. , 207203 (2014).[69] J. Knolle, D. L. Kovrizhin, J. T. Chalker, and R. Moess-ner, Phys. Rev. B , 115127 (2015).[70] J. Nasu, M. Udagawa, and Y. Motome, Phys. Rev. Lett. , 197205 (2014).[71] J. Nasu, M. Udagawa, and Y. Motome, Phys. Rev. B , 115122 (2015).[72] We average over 8 independent runs to estimate the er-rors.[73] The height of the peak in C ( T ) /L is proportional to thesystem volume, L , for first-order transitions and to L α/ν for second-order transitions.[74] M. P. M. den Nijs, Journal of Physics A: Mathematicaland General , 1857 (1979). Supplemental Material
DERIVATION OF THE MAJORANA HAMILTONIAN
The spin Hamiltonian H in Eq. (1) of the main text is exactly solvable by means of a standard procedure describedin Ref. 1 of the main text. The first step is to introduce four Majorana fermions b xj , b yj , b zj , and c j at each site j ofthe honeycomb lattice, and express the spin components σ x,y,zj in terms of these Majorana fermions as σ xj = ib xj c j , σ yj = ib yj c j , σ zj = ib zj c j . (S1)Due to the resulting enlargement of the local Hilbert space, the four Majorana fermions must be reconciled with theoriginal spin degree of freedom via the local gauge constraint − iσ xj σ yj σ zj = b xj b yj b zj c j = 1 . (S2)In turn, the corresponding gauge redundancy means that the expressions in Eq. (S1) for the spin components are notunique; for example, one can use Eq. (S2) to obtain the following equivalent expressions: σ xj = − ib yj b zj , σ yj = − ib zj b xj , σ zj = − ib xj b yj . (S3)Employing Eqs. (S1) and/or (S3) appropriately, the two terms H K and H K of the spin Hamiltonian H then become H K = − K (cid:88) (cid:104) ij (cid:105) α σ αi σ αj = − K (cid:88) (cid:104) ij (cid:105) α (cid:0) ib αi c i (cid:1)(cid:0) ib αj c j (cid:1) = iK (cid:88) (cid:104) ij (cid:105) α (cid:0) ib αi b αj (cid:1) c i c j , H K = K (cid:88) (cid:104) ijkl (cid:105) αβγ σ αi σ γj σ αk σ γl − K (cid:48) (cid:88) (cid:104) ijkl (cid:105) αβα σ αi σ γj σ γk σ αl (S4)= K (cid:88) (cid:104) ijkl (cid:105) αβγ (cid:0) ib αi c i (cid:1)(cid:0) − ib αj b βj (cid:1)(cid:0) − ib βk b γk (cid:1)(cid:0) ib γl c l (cid:1) − K (cid:48) (cid:88) (cid:104) ijkl (cid:105) αβα (cid:0) ib αi c i (cid:1)(cid:0) − ib αj b βj (cid:1)(cid:0) − ib αk b βk (cid:1)(cid:0) ib αl c l (cid:1) = iK (cid:88) (cid:104) ijkl (cid:105) αβγ (cid:0) ib αi b αj (cid:1)(cid:0) ib βk b βj (cid:1)(cid:0) ib γk b γl (cid:1) c i c l + iK (cid:48) (cid:88) (cid:104) ijkl (cid:105) αβα (cid:0) ib αi b αj (cid:1)(cid:0) ib βk b βj (cid:1)(cid:0) ib αk b αl (cid:1) c i c l . Substituting Eq. (S4) into Eq. (1) of the main text, and introducing the Z gauge fields u αij ≡ ib αi b αj , one immediatelyrecovers the Majorana Hamiltonian in Eq. (3) of the main text. Since the Z gauge fields are mutually commutingconserved quantities, u αij = ±
1, this Majorana Hamiltonian is quadratic and hence exactly solvable.
QUADRATIC MAJORANA PROBLEMS
For each phase in Fig. 2 of the main text, the ground-state flux configuration can be represented with an appropriategauge configuration u αij = ± π flux). We label eachsite of the honeycomb lattice as i ≡ (Θ , R , λ ), where Θ = { A, B } is a sublattice index, R is the lattice vector of theenlarged unit cell, and λ = 1 , . . . , n specifies the particular honeycomb unit cell within the enlarged unit cell. Notethat n = 1 for the zero-flux phase, n = 2 for the π -flux phase, and n > H = (cid:88) R , R (cid:48) (cid:88) λ,λ (cid:48) i M R (cid:48) − R ,λ,λ (cid:48) c A, R ,λ c B, R (cid:48) ,λ (cid:48) , (S5)where each M R (cid:48) − R ,λ,λ (cid:48) is a product of gauge fields u αij = ± A, R , λ ) and ( B, R (cid:48) , λ (cid:48) )occupied by the Majorana fermions c A, R ,λ and c B, R (cid:48) ,λ (cid:48) . In terms of the momentum-space complex fermions ψ A ( B ) , q ,λ = 1 √ N (cid:88) R c A ( B ) , R ,λ e − i q · R , ψ † A ( B ) , q ,λ = ψ A ( B ) , − q ,λ = 1 √ N (cid:88) R c A ( B ) , R ,λ e i q · R , (S6) a a σ v a (a) (b) (c)(d) (e)(f) (g) FIG. S1. Gauge representation of the flux configuration for the zero-flux phase (a), the π -flux phase (b), the 1 / / / / / u αij = − λ = 1 , . . . , n , each containing one A site (black dot) and one B site (white dot). where N is the number of sites, the Hamiltonian in Eq. (S5) assumes the standard Bogoliubov–de Gennes form H = 2 (cid:88) ± q (cid:88) λ,λ (cid:48) (cid:16) i ˆ M q ,λ,λ (cid:48) ψ † A, q ,λ ψ B, q ,λ (cid:48) + i ˆ M ∗ q ,λ,λ (cid:48) ψ A, q ,λ ψ † B, q ,λ (cid:48) (cid:17) = 2 (cid:88) ± q (cid:16) ψ † A, q ψ † B, q (cid:17) (cid:18) i ˆ M q − i ˆ M † q (cid:19) (cid:32) ψ A, q ψ B, q (cid:33) , ˆ M q ,λ,λ (cid:48) = (cid:88) r M r ,λ,λ (cid:48) e i q · r , ψ A ( B ) , q ≡ (cid:16) ψ A ( B ) , q , , . . . , ψ A ( B ) , q ,n (cid:17) T , (S7)where the summation is over pairs of momenta ± q due to the particle-hole redundancy ψ A ( B ) , q = ψ † A ( B ) , − q . Usingthe particular structure of this quadratic Hamiltonian, reflecting time-reversal symmetry, the Majorana energies ateach momentum ± q are then given by the singular values of the n × n matrix ˆ M q . MAJORANA NODAL STRUCTURES
From Eq. (S7), the nodal structures of the Majorana fermions are characterized by vanishing singular values of ˆ M q or, equivalently, by det ˆ M q = 0. Since det ˆ M q is generically complex, det ˆ M q = 0 translates into two independentequations for its real and imaginary parts. Consequently, for each phase in Fig. 2 of the main text, any nodal structuresare anticipated to be of codimension 2, corresponding to point nodes in two dimensions. Indeed, we generically findthat the Majorana fermions are either fully gapped or gapless at discrete points only. aR R c A , R ,1 c A , R ,2 c B , R ,1 c B , R ,2 FIG. S2. Labeling convention for the π -flux phase, including the lattice constant a , the lattice vectors R , , and the four siteswithin each unit cell. The gauge configuration is also specified; bonds with u αij = − G ( q ) as a function of the momentum components q x and q y . However, for the π -flux phase, if we only consider first-neighbor and third-neighbor Majorana hopping with ampli-tudes K and K (cid:48) , respectively, the Majorana problem takes a particularly simple form. Using the labeling conventionin Fig. S2, the matrix elements of the 2 × M q in Eq. (S7) are given byˆ M q , , = ˆ M q , , = K + K (cid:48) e − i q · ( R + R ) − K (cid:48) e i q · ( R − R ) , ˆ M q , , = K e − i q · ( R + R ) + K e − i q · R + K (cid:48) e − i q · (2 R + R ) − K (cid:48) e − i q · R + K (cid:48) + K (cid:48) e i q · ( R − R ) , (S8)ˆ M q , , = K − K e i q · R − K (cid:48) e − i q · R + K (cid:48) e i q · R + K (cid:48) e i q · ( R + R ) + K (cid:48) e i q · R , and the determinant of ˆ M q readily factorizes into the product formdet ˆ M q = K F ( q ) (cid:2) − κ (2 + κ ) − κ G ( q ) (cid:3) , (S9)where κ ≡ K (cid:48) /K is the dimensionless third-neighbor hopping amplitude, and the two functions F ( q ) and G ( q ) are F ( q ) = 1 + 2 ie − iq y a sin (cid:16) √ q x a (cid:17) ,G ( q ) = cos (cid:16) √ q x a (cid:17) − (cid:16) √ q x a (cid:17) sin (3 q y a ) . (S10)While the function F ( q ) is complex, and the solutions of F ( q ) = 0 thus give point nodes at q x a = ± π/ (6 √ πn x / √ q y a = ∓ π/ πn y / q x a = ± π/ (6 √ πn x / √ q y a = ∓ π/ πn y / n x , n y ∈ Z ), thefunction G ( q ) plotted in Fig. S3 is real, with a minimum value G min = −
3, a maximum value G max = 3 /
2, and a criticalvalue G vH = 1 corresponding to a van Hove singularity. Consequently, the solutions of 1 − κ (2 + κ ) − κ G ( q ) = 0generically correspond to nodal lines along the contours of Fig. S3 given by G ( q ) = G ( κ ) ≡ − κ (2 + κ )2 κ . (S11)To analyze these nodal lines, we plot G ( κ ) in Fig. S4 and find three critical values of κ between 0 and 1: κ = 0 . , κ = 12 (cid:0) √ − (cid:1) ≈ . , κ = 0 . . (S12)0 - / κ κ κ - κ G ( κ ) FIG. S4. Plot of the function G ( κ ) as a function of the dimensionless parameter κ = K (cid:48) /K . For κ < κ , we obtain G ( κ ) > G max , and Eq. (S11) has no solutions. For κ < κ < κ , we obtain G vH < G ( κ )
In principle, a non-magnetic impurity, such as a spin vacancy [see Fig. S5(a)], can be used as a “physical probe” todistinguish between the π -flux phases with Dirac points and Fermi surfaces of Majorana fermions. Such an impurityinduces a local modulation of the bond energy E (cid:104) ij (cid:105) α ∝ (cid:104) σ αi σ αj (cid:105) , which decays with a power law as a function ofdistance; Friedel oscillations are expected to be present (absent) in this decay if the Majorana fermions are gaplessat Fermi surfaces (Dirac points). We therefore calculate the radial Fourier transform of the bond-energy modulation,∆ E ( p ) = 1 N (cid:88) (cid:104) ij (cid:105) α e − ipr (cid:104) ij (cid:105) α (cid:104) E (cid:104) ij (cid:105) α − E (0) (cid:104) ij (cid:105) α (cid:105) , (S13)in both phases [see Fig. S5(b)], where r (cid:104) ij (cid:105) α is the distance of the bond (cid:104) ij (cid:105) α from the impurity, and E (0) (cid:104) ij (cid:105) α is thebond energy in the absence of the impurity. As expected, in the Dirac phase, ∆ E ( p ) is peaked at p = 0, while in theFermi phase, its peak is shifted to p ≈ q F , where q F is the radius of the Fermi surface [see Fig. S5(c)]. (a) -3 (b) ¯ p ¯ p ¯ p p p p κK K p p p (c) Fermi liquidDirac liquid
FIG. S5. (a) Local distortion of the Kitaev interactions ( κ (cid:54) = 1) around a nonmagnetic impurity; the particular case of a spinvacancy corresponds to κ = 0. (b) Radial Fourier transform of the bond-energy modulation around the impurity for the π -fluxphases with Dirac nodes (white circles) and Fermi surfaces (black squares) of Majorana fermions. (c) For the Fermi phase, thepeak wave vector of Friedel oscillations is compared with the characteristic dimensions of the Fermi surfaces. DYNAMICAL SPIN STRUCTURE FACTOR
The dynamical spin structure factor S µν ( q , ω ) is probed experimentally by inelastic neutron scattering. At T = 0,it is given by the spatial and temporal Fourier transform of the spin-spin correlation function in the ground state: S µν ( q , ω ) = 12 πN (cid:88) i,j (cid:90) + ∞−∞ dt e iωt − i q · ( r j − r i ) (cid:10) σ µi ( t ) σ νj (0) (cid:11) , (S14)where σ µi ( t ) ≡ e i H t σ µi e − i H t , and r i is the position of site i . For the general Hamiltonian H in Eq. (1) of the maintext, the spin-spin correlation function (cid:104) σ µi ( t ) σ νj (0) (cid:105) vanishes unless µ = ν . Moreover, its has an extremely limitedrange: (cid:104) σ µi ( t ) σ µj (0) (cid:105) is only nonzero if i and j are the same site ( i = j ) or if they are nearest-neighbor sites connectedby a µ bond (cid:104) ij (cid:105) µ . The structure factor in Eq. (S14) is thus generally given by S µν ( q , ω ) = δ µν (cid:104) S (0) µµ ( ω ) + 2 cos( q · n µ ) S (1) µµ ( ω ) (cid:105) , (S15)where S (0) µµ ( ω ) and S (1) µµ ( ω ) are on-site and nearest-neighbor contributions, S (0) µµ ( ω ) = 1 N (cid:88) i S (0) µµ,i ( ω ) , S (0) µµ,i ( ω ) = 12 π (cid:90) + ∞−∞ dt e iωt (cid:10) σ µi ( t ) σ µi (0) (cid:11) , (S16) S (1) µµ ( ω ) = 1 N (cid:88) (cid:104) ij (cid:105) µ S (1) µµ, (cid:104) ij (cid:105) µ ( ω ) , S (1) µµ, (cid:104) ij (cid:105) µ ( ω ) = 12 π (cid:90) + ∞−∞ dt e iωt (cid:10) σ µi ( t ) σ µj (0) (cid:11) , and n µ is the vector connecting the two sites i and j along any µ bond (cid:104) ij (cid:105) µ . Note that (cid:104) σ µi ( t ) σ µj (0) (cid:105) = (cid:104) σ µj ( t ) σ µi (0) (cid:105) and, consequently, S (1) µµ, (cid:104) ij (cid:105) µ ( ω ) is real due to time-reversal symmetry.The general results in Eqs. (S15) and (S16) are simplified by the unbroken space-group symmetries in each phase ofFig. 2 in the main text. First of all, due to translation symmetry, there are only a finite number of inequivalent sites i and inequivalent µ bonds (cid:104) ij (cid:105) µ , and the infinite averages in Eq. (S16) can thus be substituted with finite averagesover the n (0) inequivalent sites and the n (1) inequivalent µ bonds. In particular, n (0) = n (1) = 1 for each symmetricphase, while n (0) = n (1) = r for each flux crystal with a supercell of r plaquettes. Moreover, for all phases other thanthe 1 / z → x → y → z and therefore implies S (0 , xx ( ω ) = S (0 , yy ( ω ) = S (0 , zz ( ω ). For the 1 / y bonds, as in Fig. S1(c),the remaining point-group symmetries still imply S (0 , xx ( ω ) = S (0 , zz ( ω ).In Fig. S6, we present the dynamical spin structure factor S zz ( q , ω ) at T = 0 for each flux-crystal phase along thehigh-symmetry path M-Γ-K-M in the Brillouin zone [see inset of panel (a)]. Following the few-particle approach inRef. 69, we take the Lehmann representation of S zz ( q , ω ) and restrict our attention to intermediate states containinga single Majorana excitation. If this approximation is valid, the calculated response should approximately satisfy thesum rule (cid:82) dω S (0) zz ( ω ) = 1; for all of the results in Fig. S6, we find that (cid:82) dω S (0) zz ( ω ) > . σ zl , the Majorana fermions in the intermediate flux sector are perturbed with respect to the ground-stateflux sector, and they may even form localized states around the site l . Such a localized state corresponds to a deltapeak in the density of states and thus gives rise to a sharp feature in S zz ( q , ω ). Physically, it can be understood as amagnon bound state, σ zl = ib zl c l , of a bond Majorana fermion b zl and a matter Majorana fermion c l . MONTE CARLO IMPLEMENTATION
The Monte Carlo (MC) simulations are implemented on an L × L honeycomb lattice with dimensions L = L a and L = L a and periodic boundary conditions (PBC) in both directions (see Fig. S7). However, using these boundaryconditions, a naive numerical implementation would be extremely tedious, and we therefore simplify the procedureby removing a single bond from the lattice.Indeed, since the Hilbert space is enlarged by the introduction of the Majorana fermions, any state in the Majoranarepresentation must be projected back into the physical Hilbert space. Due to this projection, distinct states in the2 FIG. S6. Upper panel of (a)-(f): Dynamical spin structure factor S zz ( q , ω ) at T = 0 along the path M-Γ-K-M [see inset of (a)]for various phases in Fig. 2 of the main text. Each response is convolved, as a function of energy, with a Lorentzian broadeningfunction of width η = 0 . K . Lower panel of (a)-(f): Majorana density of states ρ ( ω ) for each phase. Red vertical lines are δ peaks corresponding to additional localized states that appear in the intermediate flux sector of the Lehmann representation. Majorana representation may correspond to the same physical state, and certain states in the Majorana representationmay not correspond to any physical state at all. In fact, each physical state can be represented by 2 N − distinctMajorana states, and the total fermion number, N b + N c , including bond and matter fermions, is even for all of them:( − N b ( − N c = (cid:34) (cid:89) (cid:104) jk (cid:105) α u αjk (cid:35)(cid:34) (cid:89) (cid:104) jk (cid:105) z w jk (cid:35) = +1 , (S17)where u αjk ≡ ib αj b αk (as in the main text) and w jk ≡ ic j c k . The remaining Majorana states with odd fermion numberdo not correspond to any physical state as they are annihilated by the projection. The partition function is then Z = 12 N − Tr { u αjk } even Tr { w jk } even e − β H + 12 N − Tr { u αjk } odd Tr { w jk } odd e − β H , (S18)where Tr { u αjk } even (Tr { u αjk } odd ) sums over all bond-fermion configurations { u αjk } with (cid:81) (cid:104) jk (cid:105) α u αjk = +1 ( − { w jk } even (Tr { w jk } odd ) sums over all matter-fermion configurations { w jk } with (cid:81) (cid:104) jk (cid:105) z w jk = +1 ( − (cid:104) gh (cid:105) y from the lattice (see Fig. S7) by switching off all interactions thatinvolve both sites g and h [70]. In this case, the Hamiltonian H does not depend on u ygh and, by switching between3 FIG. S7. Honeycomb lattice for the Monte Carlo simulations. We employ periodic boundary conditions in both directions butremove a single y bond (marked by the cross) to avoid a tedious numerical procedure. u ygh = ±
1, one can switch the bond-fermion parity (cid:81) (cid:104) jk (cid:105) α u αjk without affecting the spectrum of the matter fermionsat all. The partition function in Eq. (S18) can then be written as Z = 12 N − Tr { u αjk } (cid:48) (cid:104) Tr { w jk } even e − β H + Tr { w jk } odd e − β H (cid:105) = 12 N − Tr { u αjk } (cid:48) Tr { w jk } e − β H , (S19)where Tr { u αjk } (cid:48) sums over all configurations of those bond fermions that do not correspond to the removed bond (cid:104) gh (cid:105) y ,and Tr { w jk } sums over all configurations of the matter fermions. Finally, the partition function takes the form Z = 12 N − Tr { u αjk } (cid:48) (cid:89) n (cid:16) e β(cid:15) n [ { u αjk } (cid:48) ] / + e − β(cid:15) n [ { u αjk } (cid:48) ] / (cid:17) ≡ N − Tr { u αjk } (cid:48) e − βF w [ { u αjk } (cid:48) ] ≡ Z (cid:48) N − , (S20)where (cid:15) n [ { u αjk } (cid:48) ] are the non-negative single-particle energies of the matter fermions for a given configuration { u αjk } (cid:48) ofthe bond fermions. From this simplified partition function, all other thermodynamic quantities can then be obtainedas described in Ref. 70. In particular, the internal energy and the heat capacity are given by U = (cid:10) E w [ { u αjk } (cid:48) ] (cid:11) , E w [ { u αjk } (cid:48) ] ≡ − (cid:88) n (cid:15) n [ { u αjk } (cid:48) ]2 tanh β(cid:15) n [ { u αjk } (cid:48) ]2 ,C = β (cid:20)(cid:28) E w [ { u αjk } (cid:48) ] − ∂E w [ { u αjk } (cid:48) ] ∂β (cid:29) − (cid:10) E w [ { u αjk } (cid:48) ] (cid:11) (cid:21) , (S21)where (cid:104)· · · (cid:105) denotes the thermal average over all bond-fermion configurations { u αjk } (cid:48) : (cid:10) O w [ { u αjk } (cid:48) ] (cid:11) ≡ Z (cid:48) Tr { u αjk } (cid:48) (cid:104) O w [ { u αjk } (cid:48) ] e − βF w [ { u αjk } (cid:48) ] (cid:105) ..