Chiral spin liquids with crystalline ${\mathbb{Z}}_2$ gauge order in a three-dimensional Kitaev model

Chiral spin liquids (CSLs) are time-reversal-symmetry-breaking ground states of frustrated quantum magnets that show no long-range magnetic ordering but instead exhibit topological order and fractional excitations. Their realization in simple and tractable microscopic models has, however, remained an open challenge for almost two decades until it was realized that Kitaev models on lattices with odd-length loops are natural hosts for such states, even in the absence of a time-reversal-symmetry-breaking magnetic field. Here we report on the formation of CSLs in a three-dimensional Kitaev model on a hypernonagon lattice composed of nine-site loops, which differ from their widely studied two-dimensional counterparts; namely, they exhibit a crystalline ordering of the ${\mathbb{Z}}_2$ gauge fluxes and thereby break some of the underlying lattice symmetries. We study the formation of these unconventional CSLs via extensive quantum Monte Carlo simulations and demonstrate that they are separated from the featureless paramagnet at high temperatures by a single first-order phase transition at which both time-reversal and lattice symmetries are simultaneously broken. Using variational approaches for the ground state, we explore the effect of varying the Kitaev couplings and find at least five distinct CSL phases, all of which possess crystalline ordering of the ${\mathbb{Z}}_2$ gauge fluxes. For some of these phases, the complementary itinerant Majorana fermions exhibit gapless band structures with topological features such as Weyl nodes or nodal lines in the bulk and Fermi arc or drumhead surface states.

Figure 1

Schematic pictures of the Kitaev model on the hypernonagon lattice. (a) The original hypernonagon lattice [28] and (b) the distorted one for better visibility [21]. ${\mathbf{a}}_1,{\mathbf{a}}_2$, and ${\mathbf{a}}_3({\mathbf{a}}_1^{\prime },{\mathbf{a}}_2^{\prime }$, and ${\mathbf{a}}_3^{\prime })$ are the primitive translation vectors for the original (deformed) lattice. The white spheres represent the lattice sites on which the spin-$1/2$ degrees of freedom are defined. The red, green, and blue bonds denote the $x,y$, and $z$ bonds of the Kitaev model in Eq. (1), respectively. A unit cell contains 12 sites, as indicated by the numbers 1–12 in (b). Two types of inversion centers are denoted by orange and purple spheres in (b). (c)–(e) The eight types of plaquettes ${\mathrm{M}}_p(p=1$–8) in a unit cell on which the ${\mathbb{Z}}_2$ fluxes $W_p$ in Eq. (2) are defined. Each plaquette consists of nine sites.

Figure 2

3D network of the ${\mathbb{Z}}_2$ fluxes. (a) The hypernonagon lattice contains two types of closed volumes indicated by the orange and purple hexahedra. Both are constructed from the centers of the eight plaquettes, indicated by spheres. The numbers 1–8 indicate the plaquettes ${\mathrm{M}}_1$–${\mathrm{M}}_8$ in Figs. 1–1 on which the ${\mathbb{Z}}_2$ fluxes $W_p$ are defined. The spheres form a distorted cubic lattice. (b) The simple cubic lattice deformed from (a). ${\mathbf{a}}_1^{\prime },{\mathbf{a}}_2^{\prime }$, and ${\mathbf{a}}_3^{\prime }$ in (a) $[{\mathbf{a}}_1^{″},{\mathbf{a}}_2^{″}$, and ${\mathbf{a}}_3^{″}$ in (b)] are the primitive translation vectors for the deformed (simple) cubic lattice.

Figure 3

Thermodynamics of the hypernonagon Kitaev model obtained by QMC simulations for the isotropic exchange interactions $J_x=J_y=J_z=1/3$. Shown are the $T$ dependence of (a) the internal energy $\varepsilon$, (b) the specific heat $C$, and (c) the entropy $S$ per site at the last two optimization steps. The inset of (b) displays the decomposition into the contributions from the itinerant Majorana fermions and the localized ${\mathbb{Z}}_2$ fluxes. The data are calculated for the system sizes with $L=4$ and $L=6$. The gray area indicates the $T$ range plotted in Fig. 4. The dashed horizontal lines in (c) indicate 0, $\frac{1}{2}ln2$, and $ln2$.

Figure 4

Low-$T$ phase transition. Shown are the $T$ dependence of (a) the internal energy $\varepsilon$ and (b) the specific heat $C$ as given in Figs. 3 and 3, respectively, but plotted on a linear $T$ scale, zooming in the low-$T$ regime (the gray areas in Fig. 3). (c) The order parameter for the ${\mathbb{Z}}_2$ flux ordering, $S_{\mathrm{cube}}\left(\mathbf{Q}\right)$ at $\mathbf{Q}=(\pi ,\pi ,\pi )$ at the last two optimization steps. The inset of (c) indicates the flux configuration in the low-$T$ CSL phase, where black and white spheres represent fluxes with $W_p=+i$ and $-i$, respectively. ${\varepsilon }_{\mathrm{AFII}}$ in (a) indicates the ground-state energy for each system size. See the main text for details.

Figure 5

First-order nature of the thermal phase transition. Histograms of the internal energy $\varepsilon ,P\left(\varepsilon \right)$, in the vicinity of the transition temperature $T_c=0.00244\left(4\right)$ obtained by QMC simulations for (a) $L=4$ and (b) $L=6$ at the last two optimization steps. ${\varepsilon }_{\mathrm{AFII}}$ is the internal energy calculated for the $W_p$ configuration in the AFII phase for each system size. See the text for details.

Figure 6

Ground state of the hypernonagon Kitaev model. (a) The ground-state phase diagram in the parameter space $J_x+J_y+J_z=1\left(J_{\gamma }\ge 0\right)$ obtained by variational calculations. Each corner corresponds to the anisotropic limits where one of $J_{\gamma }=1$ and the other two are zero. The black dot at the center represents the isotropic case $J_x=J_y=J_z=1/3$. The white lines represent gapped/gapless boundaries with respect to the itinerant Majorana spectrum; the system has gapless excitations in the area sandwiched by the two closed white lines. Five distinct phases are found: A0F (purple), AF (magenta), AFII (gray), SI (red), and SII (yellow). The flux configurations are shown in (b)–(f) for each phase. The black and white spheres represent fluxes with $W_p=+i$ and $-i$, respectively. Since all the variational states have the same energy in the corners of the phase diagram, numerical noise may be seen close to them.

Figure 7

Gapless SI phase. (a) Evolution of Weyl nodes for coupling constants $0.267\le J_z\le 0.30$ with $J_x=J_y=(1-J_z)/2$. The sink and source indicate the negatively- and positively-charged Weyl nodes, respectively. The shaded region indicates the remaining mirror plane in the SI phase with $J_x=J_y$. (b) Corresponding Fermi arcs on the 001-surface Brillouin zone for $(J_x,J_y,J_z)=(0.36,0.36,0.28)$.

Figure 8

Gapless AF phase. (a) Evolution of Weyl nodes for coupling constants $0.42\le J_z\le 0.491$ with $J_x=J_y=(1-J_z)/2$. The sink and source indicate the negatively and positively charged Weyl nodes, respectively. (b) Corresponding Fermi arc in the 001-surface Brillouin zone for $(J_x,J_y,J_z)=(0.27,0.27,0.46)$ denoted in blue. Additional surface states are marked in red and are discussed in Appendix pp2.

Figure 9

Gapless AFII phase. (a) Nodal lines for $(J_x,J_y,J_z)=(0.34,0.34,0.32)$ shown in blue. (b) Nodal lines for $(J_x,J_y,J_z)=(0.3035,0.3035,0.393)$ shown in blue. The shaded planes denote the three mirror invariant planes, whereas the orange lines denote the $\mathbf{k}$ points invariant under a combination of mirror and inversion symmetries. Surface states for (c) $J_z=0.32$ and (d) $J_z=0.393$ depicted in blue along with projections of bulk nodal lines depicted in black for the 001-surface Brillouin zone.

Figure 10

Evolution of Weyl nodes for coupling constants $0.294\le J_z\le 0.414$ with $J_x=J_y=(1-J_z)/2$ in (a) the bulk Brillouin zone and (b) projected down to the 001-surface Brillouin zone, respectively. The two Weyl nodes on the $C_3$-invariant axis are not shown for clarity. (c) Corresponding Fermi arcs in the 001-surface Brillouin zone for $(J_x,J_y,J_z)=(0.345,0.345,0.31)$. This includes the Weyl nodes which lie on the $C_3$-invariant axis along with their corresponding Fermi arc. (d) Band structures for $J_z=0.46$. The bulk band structure along the $C_3$-invariant axis is plotted in yellow showing the two Weyl nodes. The band structure for the slab geometry along the projection of the $C_3$ axis to the surface Brillouin zone is plotted in blue and red. The projection of the two Weyl nodes as well as the Fermi arc which connects them are colored blue. The bands responsible for the remaining surface states left over from the other Weyl nodes, which are gapped out for $J_z\approx 0.414$, are colored red.