Vison Crystals in an Extended Kitaev Model on the Honeycomb Lattice

TL;DR

This paper extends the Kitaev honeycomb model with four-spin interactions, revealing a rich phase diagram with vison crystals, a $ ext{pi}$-flux spin liquid, and Lifshitz transitions, supported by Monte Carlo simulations.

Contribution

It introduces a new extended Kitaev model that maintains integrability and explores its complex phase diagram with novel vison crystal phases.

Findings

Identification of five distinct vison crystal phases.

Discovery of a symmetric $ ext{pi}$-flux spin liquid with Majorana Fermi surface.

Finite-temperature Monte Carlo results on specific heat and vison structure factor.

Abstract

We introduce an extension of the Kitaev honeycomb model by including four-spin interactions that preserve the local gauge structure and hence the integrability of the original model. The extended model has a rich phase diagram containing five distinct vison crystals, as well as a symmetric $\pi$-flux spin liquid with a Fermi surface of Majorana fermions and a sequence of Lifshitz transitions. We discuss possible experimental signatures and, in particular, present finite-temperature Monte Carlo calculations of the specific heat and the static vison structure factor. We argue that our extended model emerges naturally from generic perturbations to the Kitaev honeycomb model.

Equations

$$\mathcal{H}={\mathcal{H}}_{K_1}+{\mathcal{H}}_{K_3},$$

$${\mathcal{H}}_{K_3}=K_3\sum _{⟨ijkl{⟩}_{\alpha \beta \gamma }}{\sigma }_i^{\alpha }{\sigma }_j^{\gamma }{\sigma }_k^{\alpha }{\sigma }_l^{\gamma }-K_3^{\prime }\sum _{⟨ijkl{⟩}_{\alpha \beta \alpha }}{\sigma }_i^{\alpha }{\sigma }_j^{\gamma }{\sigma }_k^{\gamma }{\sigma }_l^{\alpha },$$

$$\mathcal{H}$$

$$+i{K}_3^{\prime }\sum _{⟨ijkl{⟩}_{\alpha \beta \alpha }}{u}_{ij}^{\alpha }{u}_{kj}^{\beta }{u}_{kl}^{\alpha }{c}_i^{}{c}_l^{},$$

$${\rho }_{\text{v}}(\mathbf{k})=\frac{1}{L^2}\sum _{p,p^{\prime }}{e}^{i\mathbf{k}\cdot ({\mathbf{X}}_p-{\mathbf{X}}_{p^{\prime }})}⟨W_p{W}_{p^{\prime }}⟩,$$

$${\sigma }_j^x=i{b}_j^x{c}_j^{},{\sigma }_j^y=i{b}_j^y{c}_j^{},{\sigma }_j^z=i{b}_j^z{c}_j^{}.$$

$$-i{\sigma }_j^x{\sigma }_j^y{\sigma }_j^z=b_j^x{b}_j^y{b}_j^z{c}_j^{}=1.$$

$${\sigma }_j^x=-i{b}_j^y{b}_j^z,{\sigma }_j^y=-i{b}_j^z{b}_j^x,{\sigma }_j^z=-i{b}_j^x{b}_j^y.$$

$${\mathcal{H}}_{K_1}$$

$${\mathcal{H}}_{K_3}$$

$$\mathcal{H}=\sum _{\mathbf{R},{\mathbf{R}}^{\prime }}\sum _{\lambda ,{\lambda }^{\prime }}i{\mathcal{M}}_{{\mathbf{R}}^{\prime }-\mathbf{R},\lambda ,{\lambda }^{\prime }}^{}{c}_{A,\mathbf{R},\lambda }^{}{c}_{B,{\mathbf{R}}^{\prime },{\lambda }^{\prime }}^{},$$

$${\psi }_{A(B),\mathbf{q},\lambda }^{}=\frac{1}{\sqrt{N}}\sum _{\mathbf{R}}{c}_{A(B),\mathbf{R},\lambda }^{}{e}^{-i\mathbf{q}\cdot \mathbf{R}},{\psi }_{A(B),\mathbf{q},\lambda }^{†}={\psi }_{A(B),-\mathbf{q},\lambda }^{}=\frac{1}{\sqrt{N}}\sum _{\mathbf{R}}{c}_{A(B),\mathbf{R},\lambda }^{}{e}^{i\mathbf{q}\cdot \mathbf{R}},$$

$$\displaystyle\mathcal{H}=2\sum_{\pm{\bm{q}}}\sum_{\lambda,\lambda^{\prime}}\left(i\hat{\mathcal{M}}_{{\bm{q}},\lambda,\lambda^{\prime}}^{\phantom{*}}\,\psi_{A,{\bm{q}},\lambda}^{{\dagger}}\,\psi_{B,{\bm{q}},\lambda^{\prime}}^{\phantom{{\dagger}}}+i\hat{\mathcal{M}}_{{\bm{q}},\lambda,\lambda^{\prime}}^{*}\,\psi_{A,{\bm{q}},\lambda}^{\phantom{{\dagger}}}\,\psi_{B,{\bm{q}},\lambda^{\prime}}^{{\dagger}}\right)=2\sum_{\pm{\bm{q}}}\left(\begin{array}[]{cc}\psi_{A,{\bm{q}}}^{{\dagger}}&\psi_{B,{\bm{q}}}^{{\dagger}}\end{array}\right)\left(\begin{array}[]{cc}0&i\hat{\mathcal{M}}_{{\bm{q}}}^{\phantom{{\dagger}}}\\ -i\hat{\mathcal{M}}_{{\bm{q}}}^{{\dagger}}&0\end{array}\right)\left(\begin{array}[]{c}\psi_{A,{\bm{q}}}^{\phantom{{\dagger}}}\\ \psi_{B,{\bm{q}}}^{\phantom{{\dagger}}}\end{array}\right),$$

$${\hat{\mathcal{M}}}_{\mathbf{q},\lambda ,{\lambda }^{\prime }}^{}=\sum _{\mathbf{r}}{\mathcal{M}}_{\mathbf{r},\lambda ,{\lambda }^{\prime }}^{}{e}^{i\mathbf{q}\cdot \mathbf{r}},{\psi }_{A(B),\mathbf{q}}^{}\equiv {\left({\psi }_{A(B),\mathbf{q},1}^{},\dots ,{\psi }_{A(B),\mathbf{q},n}^{}\right)}^T,$$

$${\hat{\mathcal{M}}}_{\mathbf{q},1,1}^{}={\hat{\mathcal{M}}}_{\mathbf{q},2,2}^{}=K_1+K_3^{\prime }{e}^{-i\mathbf{q}\cdot ({\mathbf{R}}_1+{\mathbf{R}}_2)}-K_3^{\prime }{e}^{i\mathbf{q}\cdot ({\mathbf{R}}_1-{\mathbf{R}}_2)},$$

$${\hat{\mathcal{M}}}_{\mathbf{q},1,2}^{}=K_1{e}^{-i\mathbf{q}\cdot ({\mathbf{R}}_1+{\mathbf{R}}_2)}+K_1{e}^{-i\mathbf{q}\cdot {\mathbf{R}}_2}+K_3^{\prime }{e}^{-i\mathbf{q}\cdot (2{\mathbf{R}}_1+{\mathbf{R}}_2)}-K_3^{\prime }{e}^{-i\mathbf{q}\cdot {\mathbf{R}}_1}+K_3^{\prime }+K_3^{\prime }{e}^{i\mathbf{q}\cdot ({\mathbf{R}}_1-{\mathbf{R}}_2)},$$

$${\hat{\mathcal{M}}}_{\mathbf{q},2,1}^{}=K_1-K_1{e}^{i\mathbf{q}\cdot {\mathbf{R}}_1}-K_3^{\prime }{e}^{-i\mathbf{q}\cdot {\mathbf{R}}_1}+K_3^{\prime }{e}^{i\mathbf{q}\cdot {\mathbf{R}}_2}+K_3^{\prime }{e}^{i\mathbf{q}\cdot ({\mathbf{R}}_1+{\mathbf{R}}_2)}+K_3^{\prime }{e}^{2i\mathbf{q}\cdot {\mathbf{R}}_1},$$

$$\det {\hat{\mathcal{M}}}_{\mathbf{q}}^{}=K_1^{2}F(\mathbf{q})\left[1-2\kappa (2+\kappa )-2{\kappa }^{2}G(\mathbf{q})\right],$$

$$F(\mathbf{q})$$

$$G(\mathbf{q})$$

$$G(\mathbf{q})=G_0(\kappa )\equiv \frac{1-2\kappa (2+\kappa )}{2{\kappa }^2}.$$

$$\kappa_{1}=0.2,\qquad\kappa_{2}=\frac{1}{2}\big{(}\sqrt{2}-1\big{)}\approx 0.207,\qquad\kappa_{3}=0.5.$$

$$\Delta E(p)=\frac{1}{N}\sum _{⟨ij{⟩}_{\alpha }}{e}^{-ip{r}_{⟨ij{⟩}_{\alpha }}}\left[E_{⟨ij{⟩}_{\alpha }}^{}-E_{⟨ij{⟩}_{\alpha }}^{(0)}\right],$$

$$S_{\mu\nu}({\bm{q}},\omega)=\frac{1}{2\pi N}\sum_{i,j}\int_{-\infty}^{+\infty}dt\,e^{i\omega t-i{\bm{q}}\cdot({\bm{r}}_{j}-{\bm{r}}_{i})}\,\big{\langle}\sigma_{i}^{\mu}(t)\sigma_{j}^{\nu}(0)\big{\rangle},$$

$$S_{\mu \nu }(\mathbf{q},\omega )={\delta }_{\mu \nu }\left[S_{\mu \mu }^{(0)}(\omega )+2\cos (\mathbf{q}\cdot {\mathbf{n}}_{\mu })S_{\mu \mu }^{(1)}(\omega )\right],$$

$$\displaystyle S_{\mu\mu}^{(0)}(\omega)=\frac{1}{N}\sum_{i}S_{\mu\mu,i}^{(0)}(\omega),\qquad\qquad S_{\mu\mu,i}^{(0)}(\omega)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}dt\,e^{i\omega t}\,\big{\langle}\sigma_{i}^{\mu}(t)\sigma_{i}^{\mu}(0)\big{\rangle},$$

$$\displaystyle S_{\mu\mu}^{(1)}(\omega)=\frac{1}{N}\sum_{\langle ij\rangle_{\mu}}S_{\mu\mu,\langle ij\rangle_{\mu}}^{(1)}(\omega),\qquad\,S_{\mu\mu,\langle ij\rangle_{\mu}}^{(1)}(\omega)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}dt\,e^{i\omega t}\,\big{\langle}\sigma_{i}^{\mu}(t)\sigma_{j}^{\mu}(0)\big{\rangle},$$

$$(-1)^{N_{b}}(-1)^{N_{c}}=\Bigg{[}\prod_{\langle jk\rangle_{\alpha}}u_{jk}^{\alpha}\Bigg{]}\Bigg{[}\prod_{\langle jk\rangle_{z}}w_{jk}\Bigg{]}=+1,$$

$$Z=\frac{1}{2^{N-1}}{\text{Tr}}_{\{u_{jk}^{\alpha }{\}}_{\text{even}}}{\text{Tr}}_{\{w_{jk}{\}}_{\text{even}}}{e}^{-\beta \mathcal{H}}+\frac{1}{2^{N-1}}{\text{Tr}}_{\{u_{jk}^{\alpha }{\}}_{\text{odd}}}{\text{Tr}}_{\{w_{jk}{\}}_{\text{odd}}}{e}^{-\beta \mathcal{H}},$$

$$Z=\frac{1}{2^{N-1}}\,\text{Tr}_{\{u_{jk}^{\alpha}\}^{\prime}}\Big{[}\text{Tr}_{\{w_{jk}\}_{\text{even}}}e^{-\beta\mathcal{H}}+\text{Tr}_{\{w_{jk}\}_{\text{odd}}}e^{-\beta\mathcal{H}}\Big{]}=\frac{1}{2^{N-1}}\,\text{Tr}_{\{u_{jk}^{\alpha}\}^{\prime}}\text{Tr}_{\{w_{jk}\}}e^{-\beta\mathcal{H}},$$

$$Z=\frac{1}{2^{N-1}}{\text{Tr}}_{\{u_{jk}^{\alpha }{\}}^{\prime }}\prod _n\left(e^{\beta {ϵ}_n[\{u_{jk}^{\alpha }{\}}^{\prime }]/2}+e^{-\beta {ϵ}_n[\{u_{jk}^{\alpha }{\}}^{\prime }]/2}\right)\equiv \frac{1}{2^{N-1}}{\text{Tr}}_{\{u_{jk}^{\alpha }{\}}^{\prime }}{e}^{-\beta F_w[\{u_{jk}^{\alpha }{\}}^{\prime }]}\equiv \frac{Z^{\prime }}{2^{N-1}},$$

$$U=⟨E_w[\{u_{jk}^{\alpha }{\}}^{\prime }]⟩,E_w[\{u_{jk}^{\alpha }{\}}^{\prime }]\equiv -\sum _n\frac{{ϵ}_n[\{u_{jk}^{\alpha }{\}}^{\prime }]}{2}\tanh \frac{\beta {ϵ}_n[\{u_{jk}^{\alpha }{\}}^{\prime }]}{2},$$

$$C={\beta }^2\left[⟨E_w[\{u_{jk}^{\alpha }{\}}^{\prime }{]}^2-\frac{\partial E_w[\{u_{jk}^{\alpha }{\}}^{\prime }]}{\partial \beta }⟩-{⟨E_w[\{u_{jk}^{\alpha }{\}}^{\prime }]⟩}^2\right],$$

$$\left\langle O_{w}[\{u_{jk}^{\alpha}\}^{\prime}]\right\rangle\equiv\frac{1}{Z^{\prime}}\,\text{Tr}_{\{u_{jk}^{\alpha}\}^{\prime}}\Big{[}O_{w}[\{u_{jk}^{\alpha}\}^{\prime}]\,e^{-\beta F_{w}[\{u_{jk}^{\alpha}\}^{\prime}]}\Big{]}.$$

Full text

Introduction. The famous Kitaev model on the honeycomb lattice Kitaev (2006) is an exactly solvable yet experimentally realistic model of a quantum spin liquid. In contrast to more conventional magnetic phases, quantum spin liquids retain extensive (quantum) fluctuations all the way down to zero temperature Balents (2010), where the spins appear to fractionalize into deconfined “spinon” quasiparticles coupled to appropriate gauge fields Savary and Balents (2016).

The Kitaev model is approximately realized in a family of strongly spin-orbit-coupled honeycomb materials, where its anisotropic spin interactions emerge between effective $J=1/2$ angular momenta in the $t_{2g}$ orbitals of $4d$ or $5d$ ions Jackeli and Khaliullin (2009); Rau et al. (2016); Trebst (2017); Hermanns et al. (2018). To determine the most accurate microscopic spin models for these Kitaev materials, including (Na,Li)2IrO3 Singh and Gegenwart (2010); Liu et al. (2011); Singh et al. (2012); Choi et al. (2012); Ye et al. (2012); Comin et al. (2012); Hwan Chun et al. (2015); Williams et al. (2016); Kitagawa et al. (2018) and $\alpha$-RuCl3 Plumb et al. (2014); Sandilands et al. (2015); Sears et al. (2015); Majumder et al. (2015); Johnson et al. (2015); Sandilands et al. (2016); Banerjee et al. (2016); Sears et al. (2017); Banerjee et al. (2017); Baek et al. (2017); Do et al. (2017); Banerjee et al. (2018); Hentrich et al. (2018); Kasahara et al. (2018), various extensions of the Kitaev model have been considered and analyzed with a wide range of techniques Chaloupka et al. (2010); Jiang et al. (2011); Reuther et al. (2011); Price and Perkins (2012); Rau et al. (2014); Yamaji et al. (2014); Sizyuk et al. (2014); Sela et al. (2014); Rousochatzakis et al. (2015); Kim et al. (2015); Yadav et al. (2016); Kim and Kee (2016); Winter et al. (2016); Kim et al. (2016); Janssen et al. (2016); Hou et al. (2017); Winter et al. (2017); Samarakoon et al. (2017); Ran et al. (2017); Samarakoon et al. (2018); Gordon et al. (2019). While these models are experimentally realistic and have rich phase diagrams in the classical limit, it is challenging to identify and characterize quantum phases in them. For a start, the honeycomb lattice may harbor many different quantum spin liquids Lu and Ran (2011); You et al. (2012), and the Kitaev spin liquid, captured by the Kitaev model, is only one among these many candidates. In addition, a quantum spin liquid may also remain “hidden” by appearing on top of classical symmetry-breaking order Savary and Balents (2012).

From a more phenomenological point of view, the low-energy physics of the Kitaev spin liquid is described by Majorana fermions (spinons) with Dirac nodes, coupled to an emergent $\mathbb{Z}_{2}$ gauge field Kitaev (2006). At each plaquette of the honeycomb lattice, the $\mathbb{Z}_{2}$ gauge field may form a $\pi$ flux, corresponding to a “vison” excitation. In turn, the presence of such a vison affects the kinetic energy of the spinons via the Berry phase $\pi$ picked up by each spinon moving around it. For the pure Kitaev model, the spinons are governed by a nearest-neighbor hopping problem (cf. electrons in graphene) and, due to the lack of frustration, the ground state has no visons at any plaquettes Kitaev (2006); Lieb (1994). However, if the hopping problem is frustrated by competing hopping amplitudes, the presence of a vison may reduce the frustration and thus lower the kinetic energy of the spinons. Such a frustration in the hopping amplitudes is known to stabilize crystals of topological solitons, such as baby skyrmions or merons, in itinerant magnets Ozawa et al. (2016); Batista et al. (2016); Ozawa et al. (2017), and one may thus expect it to stabilize analogous vison crystals in the Kitaev spin liquid.

In this Letter, we extend the Kitaev model by including four-spin interactions that preserve the exact solution of the model and emerge naturally from generic perturbations. By introducing frustrated further-neighbor hopping for the Majorana fermions, these additional interactions stabilize a rich variety of vison crystals, as well as a symmetric $\pi$-flux spin liquid with a vison at every plaquette. Interestingly, the $\pi$-flux spin liquid exhibits a Fermi surface of Majorana fermions undergoing two subsequent Lifshitz transitions. On a technical level, we first use a simple variational treatment to compute the zero-temperature phase diagram of our extended model. The validity of this approach is then confirmed by unbiased Monte Carlo (MC) simulations that also reveal the finite melting temperatures of the vison crystals.

Model. We consider a generalized Kitaev Hamiltonian on the honeycomb lattice:

[TABLE]

where ${\cal H}_{K_{1}}=-K_{1}\sum_{\langle ij\rangle_{\alpha}}\sigma^{\alpha}_{i}\sigma^{\alpha}_{j}$ is the usual Kitaev (2006) isotropic Kitaev Hamiltonian with ferromagnetic ($K_{1}>0$) Ising interactions between the spin components $\sigma^{\alpha}$ along each $\alpha=\{x,y,z\}$ bond $\langle ij\rangle_{\alpha}$ [see Fig. 1(a)], and

[TABLE]

where $(\alpha\beta\gamma)$ is a permutation of $(xyz)$ in each term, and $\langle ijkl\rangle_{\alpha\beta\gamma}$ is a path of length $3$ consisting of bonds $\langle ij\rangle_{\alpha}$, $\langle jk\rangle_{\beta}$, and $\langle kl\rangle_{\gamma}$. Each term in ${\cal H}_{K_{3}}$ is the product of the three terms in ${\cal H}_{K_{1}}$ that correspond to the three bonds along the appropriate path. Different $K_{3}$ and $K_{3}^{\prime}$ terms are related by space-group symmetries, simultaneously transforming the lattice and the spins; particular examples of their respective paths, with $(\alpha\beta\gamma)=(yzx)$, are depicted in Figs. 1(b) and 1(c). We remark that, for each path $\langle ijkl\rangle_{\alpha\beta\gamma}$ going around one “half” of a hexagon, connecting opposite vertices $i$ and $l$, there is a symmetry-related path $\langle lk^{\prime}j^{\prime}i\rangle_{\alpha\beta\gamma}=\langle ij^{\prime}k^{\prime}l\rangle_{\gamma\beta\alpha}$ going around the other “half” of the hexagon [see Fig. 1(b)].

Importantly, the exact solution of ${\cal H}_{K_{1}}$ Kitaev (2006) is preserved by the additional terms in Eq. (2). Indeed, since ${\cal H}$ commutes with the flux operator $W_{p}=\sigma_{1}^{x}\sigma_{2}^{y}\sigma_{3}^{z}\sigma_{4}^{x}\sigma_{5}^{y}\sigma_{6}^{z}$ at each plaquette $p$ [see Fig. 1(a)], one can identify static $\mathbb{Z}_{2}$ flux or “vison” degrees of freedom at these plaquettes, each being present (absent) if the corresponding $W_{p}$ takes eigenvalue $-1$ ($+1$). Following the Majorana fermionization $\sigma_{j}^{\alpha}=ib^{\alpha}_{j}c_{j}^{\phantom{{\dagger}}}$, the Hamiltonian takes the form sup

[TABLE]

where $u_{ij}^{\alpha}=-u_{ji}^{\alpha}\equiv ib_{i}^{\alpha}b_{j}^{\alpha}$ is a $\mathbb{Z}_{2}$ gauge field along the $\alpha$ bond $\langle ij\rangle_{\alpha}$. Since these gauge fields are conserved quantities, $u_{ij}^{\alpha}=\pm 1$, providing a redundant description of the conserved gauge fluxes, $W_{p}=u_{12}^{z}u_{32}^{x}u_{34}^{y}u_{54}^{z}u_{56}^{x}u_{16}^{y}=\pm 1$, Eq. (3) is quadratic in the Majorana fermions $c_{i}$, thus giving rise to free fermion (“spinon”) excitations after a straightforward diagonalization sup . From the perspective of the Majorana fermions, the $K_{1}$ terms describe first-neighbor hopping, while the additional $K_{3}$ and $K_{3}^{\prime}$ terms describe third-neighbor hopping not .

In analogy with how three-spin interactions may be obtained from a Zeeman field Kitaev (2006), the four-spin interactions in Eq. (2) can in principle be generated by a perturbative treatment of Heisenberg and/or symmetric off-diagonal ($\Gamma$) interactions on top of the pure Kitaev model. Taking a more universal approach and considering Eq. (3) as an effective low-energy theory for the Majorana fermions Song et al. (2016), we know that generic time-reversal-symmetric perturbations to ${\cal H}_{K_{1}}$ must generate all Majorana terms that are consistent with the projective symmetries of the Kitaev spin liquid You et al. (2012). Given that all interaction terms are irrelevant and second-neighbor hopping terms are forbidden by time reversal, Eq. (3) is the most natural effective theory beyond the pure Kitaev model.

Phase diagram. The ground state of ${\cal H}_{K_{1}}$ belongs to the zero-flux sector, characterized by $W_{p}=+1$ for all $p$ Kitaev (2006); Lieb (1994). In the presence of the additional interactions, however, the ground state may belong to a wide range of different flux sectors, as shown by the $T=0$ phase diagram in Fig. 2. This phase diagram is obtained from a simple variational analysis, by comparing the energies of the seven flux sectors appearing in the diagram on finite lattices of $48\times 48$ unit cells 111We verified that fluctuations in the phase boundaries due to finite-size effects become negligibly small for lattices larger than $36\times 36$ unit cells.. Furthermore, it is fully consistent with unbiased finite-temperature MC simulations, discussed in a later section 222We verified this statement by running unbiased MC simulations for multiple randomly chosen points within each phase on finite lattices of $12\times 12$ unit cells..

We first concentrate on the two fully symmetric non-crystal phases occupying most of the phase diagram: the zero-flux phase, which has no fluxes at any plaquettes, and the $\pi$-flux phase, which has a $\mathbb{Z}_{2}$ flux at each plaquette. For $K_{3}=K_{3}^{\prime}=0$, the creation of each $\mathbb{Z}_{2}$ flux with $W_{p}=-1$ costs a finite energy $\Delta\approx 0.15K_{1}$, and the ground state thus belongs to the zero-flux sector. For $K_{3}/K_{1}>0$, the $K_{1}$ and $K_{3}$ terms in Eq. (3) give rise to a frustrated Majorana hopping and hence an increase in the ground-state energy. However, due to the two paths between any two opposite sites $i$ and $l$ around a plaquette $p$ [see Fig. 1(b)], there are two equivalent hopping terms $\propto iK_{3}c_{i}c_{l}$ in Eq. (3), which interfere constructively for $W_{p}=+1$ and destructively for $W_{p}=-1$. Consequently, as $K_{3}/K_{1}$ is increased, fluxes are effective in relieving frustration from the Majorana hopping and thus become energetically favorable. Since the effective interaction between nearby fluxes is attractive for small $K_{3}^{\prime}/K_{1}$ Kitaev (2006), the corresponding phase transition between the zero-flux and the $\pi$-flux phases is strongly first order.

Increasing $K_{3}^{\prime}/K_{1}$, one can modify this interaction and stabilize various intermediate phases with nontrivial flux configurations. Indeed, there are five distinct translation-symmetry-breaking vison-crystal phases in Fig. 2, with their ordering wave vectors ${\bm{Q}}$ corresponding to either the $\mathrm{K}$ point or the $\mathrm{M}$ point(s) of the Brillouin zone (BZ). The two ${\bm{Q}}={\bm{Q}}_{\mathrm{K}}$ crystals have supercells of three plaquettes, containing one vison (“$1/3$ flux crystal”) and two visons (“$2/3$ flux crystal”), respectively. Since there are three different $\mathrm{M}$ points, ${\bm{Q}}={\bm{Q}}_{\mathrm{M}}$ crystals can exhibit single-${\bm{Q}}$ or multi-${\bm{Q}}$ ordering. The single-${\bm{Q}}$ crystal is a stripy configuration, corresponding to a supercell of two plaquettes containing one vison (“$1/2$ flux crystal”), while the two triple-${\bm{Q}}$ crystals have supercells of four plaquettes, containing one vison (“$1/4$ flux crystal”) and three visons (“$3/4$ flux crystal”), respectively.

Majorana problems. For the different ground-state flux sectors discussed above, distinct configurations of the gauge fields $u_{ij}^{\alpha}=\pm 1$ lead to different Majorana Hamiltonians in Eq. (3). Consequently, each phase in Fig. 2 has its own Majorana band dispersion and a corresponding density of states. The low-energy physics, giving rise to universal signatures in experiments, is determined by the nodal structures of the Majorana fermions. For the zero-flux phase, including the pure Kitaev model, as well as for the $1/4$ and $3/4$ flux crystals, the Majorana fermions are gapless at Dirac points and thus have linear density of states at low energies. For the $1/3$ and $1/2$ flux crystals, the Majorana fermions are fully gapped and thus have zero density of states below the energy gap. For the $2/3$ flux crystal, there are two disconnected phases where the Majorana fermions are gapless at Dirac points and fully gapped, respectively (see Fig. 2).

Interestingly, the Majorana fermions have more complex nodal structures in the $\pi$-flux phase. This phase is amenable to a full analytic understanding as, due to the perfect cancelation of all $K_{3}$ terms in Eq. (3), the Majorana problem has only one dimensionless parameter ratio $\kappa\equiv K_{3}^{\prime}/K_{1}$. With a simple calculation sup , we find that there are in fact three distinct $\pi$-flux phases characterized by different Majorana nodal structures.

In particular, there is a $\pi$-flux phase where the Majorana fermions are gapless at Dirac points only, and another two $\pi$-flux phases where these Dirac points coexist with Fermi surfaces (i.e., nodal lines) of distinct topologies (see Fig. 3). The dashed lines in Fig. 2 indicate two subsequent Lifshitz transitions Volovik (2003) separating these three phases as a function of increasing $\kappa$. For $\kappa<1/5$, the only nodal structures are Dirac points. At the first Lifshitz transition, $\kappa=1/5$, small pockets of Fermi surfaces appear around these Dirac points and gradually expand as $\kappa$ is further increased. At the second Lifshitz transition, $\kappa=(\sqrt{2}-1)/2\approx 0.207$, these small pockets then connect with each other to form larger pockets. We remark that the Dirac points are located at exactly the same momenta for all values of $\kappa$.

Such a coexistence of Dirac points and Fermi surfaces is rather surprising and is not expected to be stable. Instead, due to the nature of the time-reversal and particle-hole symmetries in the Majorana problem Hermanns and Trebst (2014), one would anticipate only Dirac points to be generically present, as in all the other phases of Fig. 2. Indeed, we find that the Fermi surfaces exist due to the particular simplicity of the problem up to third-neighbor hopping terms sup and that each Fermi surface is gapped out into six Dirac points (see Fig. 3) when generic fifth-neighbor hopping terms not , respecting the projective symmetries of the system, are included in Eq. (3). However, assuming that such terms are small enough, approximate Fermi surfaces are still expected to be observable in experiments.

Experimental signatures. The phase diagram in Fig. 2 contains a rich variety of phases with all possible Majorana nodal structures in two dimensions, including Fermi surfaces, Dirac points, and fully gapped scenarios. Due to their distinct low-energy physics, these phases are characterized by different experimental signatures. First, we expect the low-temperature specific heat to behave as $C\propto T$ for Fermi phases, $C\propto T^{2}$ for Dirac phases, and $C\propto e^{-\Delta_{\text{v}}/T}$ for fully gapped phases, where the activated behavior should be controlled by the vison gap $\Delta_{\text{v}}$ as it is actually smaller than the Majorana gap. Second, the various Majorana nodal structures may be distinguished by their low-energy fingerprints in spectroscopic probes, such as resonant inelastic x-ray scattering Halász et al. (2016, 2017). Third, the Majorana Fermi surface in the $\pi$-flux phase leads to impurity-induced Friedel oscillations in the magnetic energy density sup . In turn, such magnetic Friedel oscillations should be measurable with nuclear magnetic resonance (NMR) as they induce an oscillatory bond-length modulation via magnetostriction.

For the vison-crystal phases in Fig. 2, the spontaneous breaking of translation symmetry leads to further experimental signatures. First of all, due to magnetostriction, each vison crystal generates a characteristic bond-length modulation throughout the lattice, which can be picked up with NMR or elastic x-ray scattering. Moreover, the enlarged unit cell results in a larger number of distinct bands for the Majorana fermions and therefore, in contrast to the pure Kitaev model Knolle et al. (2014, 2015), the dynamical spin structure factor sup , directly measurable by inelastic neutron scattering, has multiple peaks as a function of energy (see Fig. 4). Finally, unlike the fully symmetric phases, each vison-crystal phase has a finite-temperature phase transition at a critical temperature $T_{c}$.

Monte Carlo simulations. To verify the phase diagram in Fig. 2 and to extract the melting temperatures $T_{c}$ of the vison crystals, we perform MC simulations of ${\cal H}$ based on a Metropolis algorithm to update the “classical” $\mathbb{Z}_{2}$ fields $\{u_{ij}^{\alpha}=\pm 1\}$. The energy of each field configuration is computed by diagonalizing the quadratic Majorana Hamiltonian in Eq. (3) Nasu et al. (2014, 2015) on $L\times L$ lattices with $L=\{6,12,18\}$ sup . For each temperature, a single run contains 10000 MC sweeps for equilibration and another 20000 MC sweeps for measurement 333We average over 8 independent runs to estimate the errors..

Figure 5 shows our results for the heat capacity $C(T)$ and the static vison structure factor,

[TABLE]

for representative parameters of four different vison crystals, where ${\bm{X}}_{p}$ is the position of plaquette $p$, and ${\bm{k}}$ is the ordering wave vector of each vison crystal, corresponding to either the $\mathrm{K}$ or the $\mathrm{M}$ point of the BZ. We first observe that, as for the pure Kitaev model, $C(T)$ exhibits both a high- and a low-temperature peak, which correspond to spinon and vison excitations, respectively (Nasu et al., 2014, 2015). However, the low-temperature peak signals the onset of vison-crystal ordering at $T=T_{c}$, as confirmed by the sharp growth of the corresponding Bragg peak in $\rho_{\text{v}}({\bm{k}})$. While the three lattice sizes $L=\{6,12,18\}$ do not facilitate a rigorous finite-size scaling analysis, the results in Fig. 5 suggest a first-order crystallization transition for all vison crystals, except for the $1/3$ flux crystal 444The height of the peak in $C(T)/L^{2}$ is proportional to the system volume, $L^{2}$, for first-order transitions and to $L^{\alpha/\nu}$ for second-order transitions.. Assuming that the transition into the $1/3$ flux crystal is continuous, it is conjectured to be in the universality class of the two-dimensional $3$-state Potts model, which in turn suggests that the height of the peak in $C(T)/L^{2}$ should be $\propto L^{\alpha/\nu}$ with critical exponents $\alpha=1/3$, $\nu=5/6$, and $\alpha/\nu=2/5$ den Nijs (1979). We note that, for each vison crystal, the critical temperature is $T_{c}\sim 10^{-2}K_{1}$.

Discussion. By considering a natural extension of the honeycomb Kitaev model, we have found a rich spectrum of novel spin-liquid phases that are not adiabatically connected to the original Kitaev model, including a fully symmetric $\pi$-flux spin liquid, and five distinct symmetry-breaking spin liquids with various degrees of vison crystallization. In the future, it would be interesting to study how an external magnetic field affects our spin liquids. For the Dirac phases, it may generate non-Abelian gapped spin liquids with distinct Chern numbers of the Majorana fermions Kitaev (2006). For the gapped phases, it may lead to nontrivial finite-field phase transitions between topologically distinct spin liquids.

We thank Arnab Banerjee, Hiroaki Ishizuka, and Johannes Knolle for useful comments on the manuscript. S.-S.Z., Z.W., and C.D.B. are supported by funding from the Lincoln Chair of Excellence in Physics. The work of G.B.H. at ORNL was supported by Laboratory Director’s Research and Development funds. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725.

.1 Supplemental Material

I Derivation of the Majorana Hamiltonian

The spin Hamiltonian $\mathcal{H}$ in Eq. (1) of the main text is exactly solvable by means of a standard procedure described in Ref. 1 of the main text. The first step is to introduce four Majorana fermions $b_{j}^{x}$, $b_{j}^{y}$, $b_{j}^{z}$, and $c_{j}^{\phantom{{\dagger}}}$ at each site $j$ of the honeycomb lattice, and express the spin components $\sigma_{j}^{x,y,z}$ in terms of these Majorana fermions as

[TABLE]

Due to the resulting enlargement of the local Hilbert space, the four Majorana fermions must be reconciled with the original spin degree of freedom via the local gauge constraint

[TABLE]

In turn, the corresponding gauge redundancy means that the expressions in Eq. (S1) for the spin components are not unique; for example, one can use Eq. (S2) to obtain the following equivalent expressions:

[TABLE]

Employing Eqs. (S1) and/or (S3) appropriately, the two terms ${\cal H}_{K_{1}}$ and ${\cal H}_{K_{3}}$ of the spin Hamiltonian $\mathcal{H}$ then become

[TABLE]

Substituting Eq. (I) into Eq. (1) of the main text, and introducing the $\mathbb{Z}_{2}$ gauge fields $u_{ij}^{\alpha}\equiv ib_{i}^{\alpha}b_{j}^{\alpha}$, one immediately recovers the Majorana Hamiltonian in Eq. (3) of the main text. Since the $\mathbb{Z}_{2}$ gauge fields are mutually commuting conserved quantities, $u_{ij}^{\alpha}=\pm 1$, this Majorana Hamiltonian is quadratic and hence exactly solvable.

II Quadratic Majorana problems

For each phase in Fig. 2 of the main text, the ground-state flux configuration can be represented with an appropriate gauge configuration $u_{ij}^{\alpha}=\pm 1$ (see Fig. S1). For all phases other than the zero-flux phase, the effective unit cell of the Majorana fermions is enlarged with respect to the honeycomb unit cell as a result of physical symmetry breaking (flux crystallization) and/or ostensible symmetry breaking (gauge freedom in representing each $\pi$ flux). We label each site of the honeycomb lattice as $i\equiv(\Theta,{\bm{R}},\lambda)$, where $\Theta=\{A,B\}$ is a sublattice index, ${\bm{R}}$ is the lattice vector of the enlarged unit cell, and $\lambda=1,\dots,n$ specifies the particular honeycomb unit cell within the enlarged unit cell. Note that $n=1$ for the zero-flux phase, $n=2$ for the $\pi$-flux phase, and $n>2$ for the flux-crystal phases.

Using this notation, the quadratic Majorana Hamiltonian in Eq. (3) of the main text takes the general form

[TABLE]

where each $\mathcal{M}_{{\bm{R}}^{\prime}-{\bm{R}},\lambda,\lambda^{\prime}}$ is a product of gauge fields $u_{ij}^{\alpha}=\pm 1$ along a path connecting the sites $(A,{\bm{R}},\lambda)$ and $(B,{\bm{R}}^{\prime},\lambda^{\prime})$ occupied by the Majorana fermions $c_{A,{\bm{R}},\lambda}^{\phantom{{\dagger}}}$ and $c_{B,{\bm{R}}^{\prime},\lambda^{\prime}}^{\phantom{{\dagger}}}$. In terms of the momentum-space complex fermions

[TABLE]

where $N$ is the number of sites, the Hamiltonian in Eq. (S5) assumes the standard Bogoliubov–de Gennes form

[TABLE]

where the summation is over pairs of momenta $\pm{\bm{q}}$ due to the particle-hole redundancy $\psi_{A(B),{\bm{q}}}^{\phantom{{\dagger}}}=\psi_{A(B),-{\bm{q}}}^{{\dagger}}$. Using the particular structure of this quadratic Hamiltonian, reflecting time-reversal symmetry, the Majorana energies at each momentum $\pm{\bm{q}}$ are then given by the singular values of the $n\times n$ matrix $\hat{\mathcal{M}}_{{\bm{q}}}^{\phantom{{\dagger}}}$.

III Majorana nodal structures

From Eq. (S13), the nodal structures of the Majorana fermions are characterized by vanishing singular values of $\hat{\mathcal{M}}_{{\bm{q}}}^{\phantom{{\dagger}}}$ or, equivalently, by $\det\hat{\mathcal{M}}_{{\bm{q}}}^{\phantom{{\dagger}}}=0$. Since $\det\hat{\mathcal{M}}_{{\bm{q}}}^{\phantom{{\dagger}}}$ is generically complex, $\det\hat{\mathcal{M}}_{{\bm{q}}}^{\phantom{{\dagger}}}=0$ translates into two independent equations for its real and imaginary parts. Consequently, for each phase in Fig. 2 of the main text, any nodal structures are anticipated to be of codimension $2$, corresponding to point nodes in two dimensions. Indeed, we generically find that the Majorana fermions are either fully gapped or gapless at discrete points only.

However, for the $\pi$-flux phase, if we only consider first-neighbor and third-neighbor Majorana hopping with amplitudes $K_{1}$ and $K_{3}^{\prime}$, respectively, the Majorana problem takes a particularly simple form. Using the labeling convention in Fig. S2, the matrix elements of the $2\times 2$ matrix $\hat{\mathcal{M}}_{{\bm{q}}}^{\phantom{{\dagger}}}$ in Eq. (S13) are given by

[TABLE]

and the determinant of $\hat{\mathcal{M}}_{{\bm{q}}}^{\phantom{{\dagger}}}$ readily factorizes into the product form

[TABLE]

where $\kappa\equiv K_{3}^{\prime}/K_{1}$ is the dimensionless third-neighbor hopping amplitude, and the two functions $F({\bm{q}})$ and $G({\bm{q}})$ are

[TABLE]

While the function $F({\bm{q}})$ is complex, and the solutions of $F({\bm{q}})=0$ thus give point nodes at $q_{x}a=\pm\pi/(6\sqrt{3})+2\pi n_{x}/\sqrt{3}$ and $q_{y}a=\mp\pi/6+2\pi n_{y}/3$ as well as at $q_{x}a=\pm 5\pi/(6\sqrt{3})+2\pi n_{x}/\sqrt{3}$ and $q_{y}a=\mp\pi/6+2\pi n_{y}/3$ (with $n_{x},n_{y}\in\mathbb{Z}$), the function $G({\bm{q}})$ plotted in Fig. S3 is real, with a minimum value $G_{\mathrm{min}}=-3$, a maximum value $G_{\mathrm{max}}=3/2$, and a critical value $G_{\mathrm{vH}}=1$ corresponding to a van Hove singularity. Consequently, the solutions of $1-2\kappa(2+\kappa)-2\kappa^{2}G({\bm{q}})=0$ generically correspond to nodal lines along the contours of Fig. S3 given by

[TABLE]

To analyze these nodal lines, we plot $G_{0}(\kappa)$ in Fig. S4 and find three critical values of $\kappa$ between [math] and $1$:

[TABLE]

For $\kappa<\kappa_{1}$, we obtain $G_{0}(\kappa)>G_{\mathrm{max}}$, and Eq. (S17) has no solutions. For $\kappa_{1}<\kappa<\kappa_{2}$, we obtain $G_{\mathrm{vH}}<G_{0}(\kappa)<G_{\mathrm{max}}$, and the solutions of Eq. (S17) are nodal lines surrounding the maxima of $G({\bm{q}})$. Interestingly, these maxima coincide with the point nodes characterized by $F({\bm{q}})=0$. Finally, for $\kappa>\kappa_{2}$, we obtain $G_{\mathrm{min}}\leq G_{0}(\kappa)<G_{\mathrm{vH}}$, and the solutions of Eq. (S17) are nodal lines surrounding the minima of $G({\bm{q}})$. In particular, for $\kappa=\kappa_{3}$, these nodal lines contract to point nodes as $G_{0}(\kappa)=G_{\mathrm{min}}$.

IV Magnetic Friedel oscillations

In principle, a non-magnetic impurity, such as a spin vacancy [see Fig. S5(a)], can be used as a “physical probe” to distinguish between the $\pi$-flux phases with Dirac points and Fermi surfaces of Majorana fermions. Such an impurity induces a local modulation of the bond energy $E_{\langle ij\rangle_{\alpha}}\propto\langle\sigma_{i}^{\alpha}\sigma_{j}^{\alpha}\rangle$, which decays with a power law as a function of distance; Friedel oscillations are expected to be present (absent) in this decay if the Majorana fermions are gapless at Fermi surfaces (Dirac points). We therefore calculate the radial Fourier transform of the bond-energy modulation,

[TABLE]

in both phases [see Fig. S5(b)], where $r_{\langle ij\rangle_{\alpha}}$ is the distance of the bond $\langle ij\rangle_{\alpha}$ from the impurity, and $E_{\langle ij\rangle_{\alpha}}^{(0)}$ is the bond energy in the absence of the impurity. As expected, in the Dirac phase, $\Delta E(p)$ is peaked at $p=0$, while in the Fermi phase, its peak is shifted to $p\approx 2q_{F}$, where $q_{F}$ is the radius of the Fermi surface [see Fig. S5(c)].

V Dynamical spin structure factor

The dynamical spin structure factor $S_{\mu\nu}({\bm{q}},\omega)$ 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:

[TABLE]

where $\sigma_{i}^{\mu}(t)\equiv e^{i\mathcal{H}t}\sigma_{i}^{\mu}e^{-i\mathcal{H}t}$, and ${\bm{r}}_{i}$ is the position of site $i$. For the general Hamiltonian $\mathcal{H}$ in Eq. (1) of the main text, the spin-spin correlation function $\langle\sigma_{i}^{\mu}(t)\sigma_{j}^{\nu}(0)\rangle$ vanishes unless $\mu=\nu$. Moreover, its has an extremely limited range: $\langle\sigma_{i}^{\mu}(t)\sigma_{j}^{\mu}(0)\rangle$ is only nonzero if $i$ and $j$ are the same site ($i=j$) or if they are nearest-neighbor sites connected by a $\mu$ bond $\langle ij\rangle_{\mu}$. The structure factor in Eq. (S20) is thus generally given by

[TABLE]

where $S_{\mu\mu}^{(0)}(\omega)$ and $S_{\mu\mu}^{(1)}(\omega)$ are on-site and nearest-neighbor contributions,

[TABLE]

and ${\bm{n}}_{\mu}$ is the vector connecting the two sites $i$ and $j$ along any $\mu$ bond $\langle ij\rangle_{\mu}$. Note that $\langle\sigma_{i}^{\mu}(t)\sigma_{j}^{\mu}(0)\rangle=\langle\sigma_{j}^{\mu}(t)\sigma_{i}^{\mu}(0)\rangle$ and, consequently, $S_{\mu\mu,\langle ij\rangle_{\mu}}^{(1)}(\omega)$ is real due to time-reversal symmetry.

The general results in Eqs. (S21) and (S22) are simplified by the unbroken space-group symmetries in each phase of Fig. 2 in the main text. First of all, due to translation symmetry, there are only a finite number of inequivalent sites $i$ and inequivalent $\mu$ bonds $\langle ij\rangle_{\mu}$, and the infinite averages in Eq. (S22) can thus be substituted with finite averages over the $n_{(0)}$ inequivalent sites and the $n_{(1)}$ inequivalent $\mu$ bonds. In particular, $n_{(0)}=n_{(1)}=1$ for each symmetric phase, while $n_{(0)}=n_{(1)}=r$ for each flux crystal with a supercell of $r$ plaquettes. Moreover, for all phases other than the $1/2$ flux crystal, there is a threefold rotation symmetry around the center of some plaquette, which permutes the spin components as $z\rightarrow x\rightarrow y\rightarrow z$ and therefore implies $S_{xx}^{(0,1)}(\omega)=S_{yy}^{(0,1)}(\omega)=S_{zz}^{(0,1)}(\omega)$. For the $1/2$ flux crystal, this rotation symmetry is spontaneously broken; if the flux stripes are perpendicular to the $y$ bonds, as in Fig. S1(c), the remaining point-group symmetries still imply $S_{xx}^{(0,1)}(\omega)=S_{zz}^{(0,1)}(\omega)$.

In Fig. S6, we present the dynamical spin structure factor $S_{zz}({\bm{q}},\omega)$ at $T=0$ for each flux-crystal phase along the high-symmetry path $\mathrm{M}$-$\Gamma$-$\mathrm{K}$-$\mathrm{M}$ in the Brillouin zone [see inset of panel (a)]. Following the few-particle approach in Ref. Knolle et al., 2015, we take the Lehmann representation of $S_{zz}(\bm{q},\omega)$ and restrict our attention to intermediate states containing a single Majorana excitation. If this approximation is valid, the calculated response should approximately satisfy the sum rule $\int d\omega\,S_{zz}^{(0)}(\omega)=1$; for all of the results in Fig. S6, we find that $\int d\omega\,S_{zz}^{(0)}(\omega)>0.6$.

The energy dependence of the dynamical spin structure factor reflects the Majorana density of states in the intermediate flux sector of the Lehmann representation (see Fig. S6). Due to the two fluxes created (or destroyed) by the spin operator $\sigma_{l}^{z}$, the Majorana fermions in the intermediate flux sector are perturbed with respect to the ground-state flux sector, and they may even form localized states around the site $l$. Such a localized state corresponds to a delta peak in the density of states and thus gives rise to a sharp feature in $S_{zz}({\bm{q}},\omega)$. Physically, it can be understood as a magnon bound state, $\sigma_{l}^{z}=ib_{l}^{z}c_{l}^{\phantom{{\dagger}}}$, of a bond Majorana fermion $b_{l}^{z}$ and a matter Majorana fermion $c_{l}^{\phantom{{\dagger}}}$.

VI Monte Carlo implementation

The Monte Carlo (MC) simulations are implemented on an $L\times L$ honeycomb lattice with dimensions $\bm{L}_{1}=L\bm{a}_{1}$ and $\bm{L}_{2}=L\bm{a}_{2}$ and periodic boundary conditions (PBC) in both directions (see Fig. S7). However, using these boundary conditions, a naive numerical implementation would be extremely tedious, and we therefore simplify the procedure by 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 Majorana representation must be projected back into the physical Hilbert space. Due to this projection, distinct states in the Majorana representation may correspond to the same physical state, and certain states in the Majorana representation may not correspond to any physical state at all. In fact, each physical state can be represented by $2^{N-1}$ distinct Majorana states, and the total fermion number, $N_{b}+N_{c}$, including bond and matter fermions, is even for all of them:

[TABLE]

where $u_{jk}^{\alpha}\equiv ib_{j}^{\alpha}b_{k}^{\alpha}$ (as in the main text) and $w_{jk}\equiv ic_{j}c_{k}$. The remaining Majorana states with odd fermion number do not correspond to any physical state as they are annihilated by the projection. The partition function is then

[TABLE]

where $\text{Tr}_{\{u_{jk}^{\alpha}\}_{\text{even}}}$ ($\text{Tr}_{\{u_{jk}^{\alpha}\}_{\text{odd}}}$) sums over all bond-fermion configurations $\{u_{jk}^{\alpha}\}$ with $\prod_{\langle jk\rangle_{\alpha}}u_{jk}^{\alpha}=+1$ ($-1$), while $\text{Tr}_{\{w_{jk}\}_{\text{even}}}$ ($\text{Tr}_{\{w_{jk}\}_{\text{odd}}}$) sums over all matter-fermion configurations $\{w_{jk}\}$ with $\prod_{\langle jk\rangle_{z}}w_{jk}=+1$ ($-1$).

For the relatively small system sizes accessible with MC, it is important to take proper care of the even/odd projections, which in turn makes the MC implementation extremely tedious. However, the computation is simplified tremendously if we remove a single bond $\langle gh\rangle_{y}$ from the lattice (see Fig. S7) by switching off all interactions that involve both sites $g$ and $h$ Nasu et al. (2014). In this case, the Hamiltonian $\mathcal{H}$ does not depend on $u_{gh}^{y}$ and, by switching between $u_{gh}^{y}=\pm 1$, one can switch the bond-fermion parity $\prod_{\langle jk\rangle_{\alpha}}u_{jk}^{\alpha}$ without affecting the spectrum of the matter fermions at all. The partition function in Eq. (S24) can then be written as

[TABLE]

where $\text{Tr}_{\{u_{jk}^{\alpha}\}^{\prime}}$ sums over all configurations of those bond fermions that do not correspond to the removed bond $\langle gh\rangle_{y}$, and $\text{Tr}_{\{w_{jk}\}}$ sums over all configurations of the matter fermions. Finally, the partition function takes the form

[TABLE]

where $\epsilon_{n}[\{u_{jk}^{\alpha}\}^{\prime}]$ are the non-negative single-particle energies of the matter fermions for a given configuration $\{u_{jk}^{\alpha}\}^{\prime}$ of the bond fermions. From this simplified partition function, all other thermodynamic quantities can then be obtained as described in Ref. Nasu et al., 2014. In particular, the internal energy and the heat capacity are given by

[TABLE]

where $\langle\cdots\rangle$ denotes the thermal average over all bond-fermion configurations $\{u_{jk}^{\alpha}\}^{\prime}$:

[TABLE]

Note that the same simplification in the partition function could also be achieved by removing several bonds from the lattice; we remove only one bond to minimize any accompanying finite-size effects.