Fracton topological order from nearest-neighbor two-spin interactions and dualities

Fracton topological order describes a remarkable phase of matter, which can be characterized by fracton excitations with constrained dynamics and a ground-state degeneracy that increases exponentially with the length of the system on a three-dimensional torus. However, previous models exhibiting this order require many-spin interactions, which may be very difficult to realize in a real material or cold atom system. In this work, we present a more physically realistic model which has the so-called X-cube fracton topological order [Vijay, Haah, and Fu, Phys. Rev. B 94, 235157 (2016)] but only requires nearest-neighbor two-spin interactions. The model lives on a three-dimensional honeycomb-based lattice with one to two spin-1/2 degrees of freedom on each site and a unit cell of six sites. The model is constructed from two orthogonal stacks of $Z_2$ topologically ordered Kitaev honeycomb layers [Kitaev, Ann. Phys.321, 2 (2006)], which are coupled together by a two-spin interaction. It is also shown that a four-spin interaction can be included to instead stabilize 3+1D $Z_2$ topological order. We also find dual descriptions of four quantum phase transitions in our model, all of which appear to be discontinuous first-order transitions.

Figure 1

Two series of condensation driven phase transitions with dual descriptions [(b) and (c)], which are described in Sec. 3 and Sec. 4c. All four phase transitions appear to be first order when generic perturbations are added (Appendix pp3). The transition from XY and XZ plane 2+1D $Z_2$ quantum spin liquid (QSL) to X-cube fracton topological order is driven by a condensation of YZ-plane composite flux loops [red loop in (b)], which are excitations of a loop of XY and XZ–plane 2+1D toric code flux operators $({\left[\tilde{\sigma }\right]}_p^{{\square }_{\text{XY}}}$ and ${\left[\tilde{\tau }\right]}_p^{{\square }_{\text{XZ}}}$ in Sec. 3a) [21]. In the dual theory, the YZ-plane flux loop is mapped to a YZ-plane domain wall excitation of Ising spins [green cones in (b), ${\eta }_{\widehat{\iota }}^z$ in Sec. 3a]. The transition from X-cube order to a YZ-plane 2+1D $Z_2$ QSL is driven by a condensation of dimension-1 X-cube particle excitations [red cross in (c), ${\left[\widehat{\rho }\right]}_{\iota }^{{+}_{\text{XY}}}$ and ${\left[\widehat{\rho }\right]}_{\iota }^{{+}_{\text{XZ}}}$ in Sec. 3b], which are bound to 1D lines parallel to the $x$ axis. In the dual theory, these $x$-axis particle excitations are mapped to excitations [light green cone $|\leftarrow \rangle$ in (c), ${\mu }_{\iota }^x$ in Sec. 3b] of 1+1D Ising paramagnets, which spontaneously break a $Z_2$ Ising symmetry when they condense since they are conserved modulo two. The transition from XY & XZ plane QSL to 3+1D QSL is driven by condensing a composite of XY and XZ–plane 2+1D toric code charge excitations [red cross in (c)], which is also dual to a 1+1D Ising paramagnet excitation. The transition from 3+1D QSL to YZ-plane QSL is driven by condensing YZ-plane flux loops [red in (b)] in 3+1D toric code, which are dual to 2+1D Ising domain walls. These dualities also describe phase transitions shown in Figs. 3 and 8.

Figure 2

The lattice that we consider: orthogonal stacks of honeycomb lattices along the XY and XZ planes which intersect on the dark-blue links, which we refer to as shared links in the text. ${\sigma }^{\mu }$ and ${\tau }^{\mu }$ spin-1/2 degrees of freedom reside on sites of the XY and XZ honeycomb lattices, respectively. Thus the shared sites on either end of a dark-blue link host both a ${\sigma }^{\mu }$ and ${\tau }^{\mu }$ spin. The unit cell consists of eight spins: the spins on both sides of a $\sigma$ link, a $\tau$ link, and a shared link, which has two spins on each side. On each honeycomb lattice, we impose a Kitaev honeycomb model [Eq. (2)] where the links are colored red $(x$ link), green $(y$ link), or (dark) blue $(z$ link) to indicate the kind of coupling in Eq. (2): ${\sigma }_i^x{\sigma }_j^x,{\sigma }_i^y{\sigma }_j^y,{\sigma }_i^z{\sigma }_j^z$, respectively (and similar for ${\tau }^{\mu })$. To induce X-cube fracton topological order, we place a strong ${\sigma }_i^z{\tau }_i^z$ coupling on every shared site.

Figure 3

Rough phase diagram of the Hamiltonian in Eq. (1) as a function of $t$ (horizontal axis) and $K_x=K_y$ (vertical axis); $K_z$ will be used to set the energy scale. When $t=0$, the model consists of two orthogonal stacks (along the XY and XZ planes) of decoupled honeycomb lattices (Fig. 2), where each honeycomb lattice will host a quantum spin liquid (QSL) with $Z_2$ topological order, which will be gapped when $K_x=K_y<K_z/2$ and gapless when $K_x=K_y\ge K_z/2$ [20]. These QSL phases are stable to a small coupling $t$ which couples the orthogonal honeycomb lattices together. (The gapless QSL is stable since $t$ preserves time reversal symmetry; see, e.g., Sec. 6.1 of Ref. [20].) In Sec. 2b and Appendix pp1-s2, we argue that X-cube fracton topological order [16] results when $max(K_x,K_y)\ll min(K_z,t)$. We do not know what phase(s) occur in the white region. In the bottom left corner, a phase transition between the QSL and fracton orders occurs near $t\sim K_x^2{K}_y^2/K_z^3$ [in the $max(K_x,K_y)\ll K_z$ limit], which can be inferred from Fig. 11 in the appendix. In Sec. 3a, we propose that this transition is dual to a stack 2+1D Ising transitions; in Appendix pp3, we find that the transition is likely first order and that (depending on the sign of certain perturbations) a topologically ordered intermediate phase is also possible.

Figure 4

An XY plane honeycomb lattice: a constant $z$ slice of Fig. 2. There is a ${\sigma }^{\mu }$ degree of freedom on every site, and also a ${\tau }^{\mu }$ on every shared site: i.e., every site next to a dark-blue shared link, which overlap gray lines. Red, green, and blue links host a $K_x{\sigma }_i^x{\sigma }_j^x,K_y{\sigma }_i^y{\sigma }_j^y$, or $K_z{\sigma }_i^z{\sigma }_j^z$ coupling, respectively. When $max(K_x,K_y)\ll K_z$, degenerate perturbation theory will produce the operators ${\left[\sigma \right]}_p^{{\square }_{\text{XY}}}$ (top left) and ${\left[\sigma \right]}_{\iota }^{{+}_{\text{XY}}}$ (top right), which both take the form ${\sigma }_{i_1}^z{\sigma }_{i_2}^y{\sigma }_{i_3}^x{\sigma }_{i_4}^z{\sigma }_{i_5}^y{\sigma }_{i_6}^x$ for sites around a hexagon, as shown above. The low-energy Hilbert space will have ${\sigma }_i^z{\sigma }_j^z=1$ across all blue links; and so we will perform a local change of basis (and projection) $\sigma \to \tilde{\sigma }$ given in Eq. (5) to describe the effective low-energy Hilbert space; the new ${\tilde{\sigma }}^{\mu }$ operators are centered on the blue links. In the $\tilde{\sigma }$ basis, ${\left[\sigma \right]}_p^{{\square }_{\text{XY}}}\to {\left[\tilde{\sigma }\right]}_p^{{\square }_{\text{XY}}}$ (bottom left) and ${\left[\sigma \right]}_{\iota }^{{+}_{\text{XY}}}\to {\left[\tilde{\sigma }\right]}_{\iota }^{{+}_{\text{XY}}}$ (bottom right) take the form of toric code flux and charge operators, as shown above. The gray lines indicate the rectangular lattice of the toric code. The XZ planes are similar, but with the replacement $\sigma \leftrightarrow \tau$. See Fig. 5 for a three-dimensional version.

Figure 5

The lattice that we consider (Fig. 2), except we also draw gray lines (as in Fig. 4) to indicate the intersecting rectangular lattices of the toric code [Eq. (4)] on XY and XZ planes.

Figure 6

The $Z_2$ fracton number operator generated from degenerate perturbation theory [Eq. (12)] when $max(K_x,K_y)\ll min(K_z,t)$. The operator is expressed in three different basis: is in original basis [Eq. (1)], is in the 2D toric code basis [Eq. (5)], and is in the X-cube model basis [Eq. (14)]. is the product of two XY-plane and two XZ-plane plaquette/flux operators $({\left[\tilde{\sigma }\right]}_p^{{\square }_{\text{XY}}}$ and ${\left[\tilde{\tau }\right]}_p^{{\square }_{\text{XZ}}}$, Fig. 4) around the cube.

Figure 7

The ${\left[\widehat{\rho }\right]}_{\iota }^{{+}_{\text{XY}}},{\left[\widehat{\rho }\right]}_{\iota }^{{+}_{\text{XZ}}}$, and ${\left[\widehat{\rho }\right]}_{\iota }^{{+}_{\text{YZ}}}$ operators in the X-cube basis [Eq. (14)]. These operators, along with (Fig. 6) are the commuting operators used to define the X-cube fracton Hamiltonian [Eq. (12) after a change of basis (14)]. ${\left[\sigma \right]}_{\iota }^{{+}_{\text{XY}}}$ was given in the original basis in Fig. 4; ${\left[\tau \right]}_{\iota }^{{+}_{\text{XZ}}}$ is similar, but with $\sigma \leftrightarrow \tau$ and on a XZ plane instead of a XY plane. ${[\sigma ,\tau ]}_{\iota }^{{+}_{\text{YZ}}}$ is given by ${[\sigma ,\tau ]}_{\iota }^{{+}_{\text{YZ}}}={\left[\sigma \right]}_{\iota }^{{+}_{\text{XY}}}{\left[\tau \right]}_{\iota }^{{+}_{\text{XZ}}}$.

Figure 8

Rough phase diagram of $H^{\prime }$ [Eq. (25)] as a function of $t$ (horizontal axis) and $h$ (vertical axis) in the limit $max(K_x,K_y)\ll K_z$ with $max(K_x,K_y)/K_z$ small but finite and $h^{\prime }$ infinitesimally small (to break accidental degeneracies). In Fig. 11, we utilize infinitesimal $max(K_x,K_y)/K_z$ in order display more precise claims regarding phase boundaries. Rough phase boundaries are computed in Appendix pp1. Four phase transitions occur around the edge of the phase diagram. In Secs. 3 and 4c, we present dual descriptions of these four (likely first-order) phase transitions. We do not know what phase(s) occur in the white region; this “region of ignorance” grows as we move away from the $max(K_x,K_y)\ll K_z$ limit. In Appendix pp3, we find that (depending on the sign of certain perturbations) topologically ordered intermediate phases are also possible in between all of the above phase transitions.

Figure 9

The YZ-plane toric code plaquette/flux operator expressed in three different basis: ${[\sigma ,\tau ]}_p^{{\square }_{\text{YZ}}}$ (top) is in the original basis [Eq. (1)], ${[\tilde{\sigma },\tilde{\tau }]}_p^{{\square }_{\text{YZ}}}$ (mid-left) is in the 2D toric code basis [Eq. (5)], and ${\left[\widehat{\rho }\right]}_p^{{\square }_{\text{YZ}}}$ (mid-right) is the X-cube basis [Eq. (14)]. ${\left[\widehat{\rho }\right]}_{\iota }^{*}$ (bottom) is the vertex/charge operator in 3+1D toric code [Eq. (33)].

Figure 10

Phase diagram of ${\tilde{H}}_{\text{global}}$ [Eq. (43)] for $E_{\square }^{\text{YZ}}\ll E_{+}$. ${\tilde{H}}_{\text{global}}$ is used to derive dual descriptions of the phase transitions between the above phases, which occur on the lines $t\sim E_{\square }$ and $h=2E_{+}$. The phase boundaries in Fig. 8 are different because that phase diagram describes a model $H^{\prime }$ [Eq. (25)] where the effective $E_{\square }$ and $E_{+}$ depend on $t$ and $h$ as shown in Eq. (A1).

Figure 11

Phase diagram of $H^{\prime }$ [Eq. (25)] in the limit $max(K_x,K_y)\ll K_z=1$ with $max(K_x,K_y)/K_z$ infinitesimally small. Finite $max(K_x,K_y)/K_z$ is shown in Figs. 3 and 8. $Z_2$ QSL is short for quantum spin liquid with $Z_2$ topological order. The 2+1D QSL phases are in the same phase as parallel layers of 2+1D toric code along certain orthogonal planes. The phase boundaries are calculated in Appendix pp1 and given by Eqs. (A2, A3, A4, A5). Degenerate perturbation theory does not determine the nature of the phase transitions, nor the possible existence of intermediate phases. However, duality arguments can be used to describe these transitions (Secs. 3 and 4c). The dotted gray lines separate overlapping regions where different choices of unperturbed Hamiltonian $H_0$ were used in degenerate perturbation theory in Appendix pp1. The diagram is a log-log plot where the axes have log base $K_{xy}^{-1}$ where $K_{xy}=\sqrt{K_x{K}_y}$. (The axes are labeled by powers of $K_{xy}$, which is to be thought of as a very small number $\epsilon$.) Note that the point set topology of the log-log coordinates in the $K_{xy}\ll K_z$ limit can be surprising. Specifically, for finite $K_{xy}/K_z$, intermediate phases are possible (see Appendix pp3), but in the $K_{xy}\ll K_z$ limit, these intermediate phases shrink into the red lines.

Figure 12

The first column of plots shows ${\xi }_c/L$ vs the ferromagnetic Ising coupling $J$, where ${\xi }_c$ is the connected correlation length between two spins on the same layer. The second column is a closeup of a phase transition. The third column is the squared layer magnetization $\langle S_L^2\rangle$ vs $J$, where $S_L\equiv \frac{1}{N_I}{\sum }_I{\sigma }_{LI}$ is the average magnetization on a given layer. The system size is $L\times L\times L\times L$. Each row corresponds to a different dimension for each layer, or different sign of $K$. All phase transitions appear to be first order since at each transition, ${\xi }_c/L$ tends to decrease as the system length $L$ increases. An exception is the second transition for the 2D stack of 1+1D layers with $K=-1/8$ (fourth row and firdt column of plots, near $J\sim 0.88)$. However, when $K$ is increased to $K=-1/4$, the transition looks first order. (When $K=-1/4$, we cannot accurately measure the correlation length below $J\sim 1.2$ due to large autocorrelation times in this region.) We expect that larger system sizes would show that the $K=-1/8$ transition is also second order. Some of the data is a little messy (e.g., not smooth vs $J)$ due to a loss of ergodicity (due to large autocorrelation times) near the first-order transitions, which is a well known shortcoming of the Wolff MC update algorithm near a first-order phase transition. (The one standard deviation statistical error bars do not account for this error.) All data points used at least $2^{12}$ MC sweeps. The second row used $2^{16}$ sweeps.

Figure 13

We show that the order parameter $R$ [Eq. (C3)] describes the intermediate phase.