Fracton Topological Order from Nearest-Neighbor Two-Spin Interactions and Dualities

TL;DR

This paper introduces a realistic 3D lattice model with only nearest-neighbor two-spin interactions that exhibits fracton topological order, making it more feasible for experimental realization.

Contribution

The authors present a new 3D lattice model with only two-spin interactions that realizes X-cube fracton order, derived from coupled Kitaev honeycomb layers, and analyze its phase transitions.

Findings

Model exhibits X-cube fracton topological order with simple interactions.

Includes dual descriptions of quantum phase transitions.

Demonstrates first-order nature of phase transitions.

Abstract

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 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 6 sites. The model is constructed from two orthogonal stacks of $Z_2$ topologically ordered Kitaev honeycomb layers, 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.

Figures

Tables

Table 1. Table 1: Ground state degeneracies for the four phases in Fig.  1 , along with the four possible intermediate phases that can occur. For example, intermediate #3 can occur between the 2+1D Z 2 subscript 𝑍 2 Z_{2} QSL and X-cube phases. The degeneracies of all of the intermediate phases take the form 2 × ⋯ 2 ⋯ 2\times\cdots . The first factor of two in the degeneracies is not topological, and occurs due to the spontaneous translation symmetry breaking of the intermediate phase. The rest of the degeneracy is topological.

	XY&XZ plane				YZ-plane
phase	2+1D $Z_2$ QSL	intermediate #1	3+1D $Z_2$ QSL	intermediate #2	2+1D $Z_2$ QSL
degeneracy	$2^{2L_y+2L_z}$	$2\times 2^6$	$2^3$	$2\times 2^{L_x}$	$2^{2L_x}$
phase		intermediate #3	X-cube	intermediate #4
degeneracy		$2\times 2^{L_x+2L_y+2L_z-3}$	$2^{2L_x+2L_y+2L_z-3}$	$2\times 2^{2L_x+3}$

Equations

$$H$$

$$H_K$$

$$H_t$$

$$H_{\text{toric}}$$

$$H_{\text{toric}}^{\text{XY}}$$

$$H_{\text{toric}}^{\text{XZ}}$$

$$E_{\square }\sim E_{+}$$

$$[\sigma {]}_p^{{\square }_{\text{XY}}},[\tau {]}_p^{{\square }_{\text{XZ}}},[\sigma {]}_{\iota }^{{+}_{\text{XY}}},[\tau {]}_{\iota }^{{+}_{\text{XZ}}}\text{are defined in Fig.reffig:toricCode}$$

$${\sigma }_i^z{\sigma }_j^z$$

$${\sigma }_i^z{\sigma }_j^0\&{\sigma }_i^0{\sigma }_j^z$$

$${\sigma }_i^x{\sigma }_j^y$$

$$\tau$$

$$H_{\text{toric}}$$

$${\tilde{H}}_{\text{toric}}^{\text{XY}}$$

$${\tilde{H}}_{\text{toric}}^{\text{XZ}}$$

$$[\sigma {]}_p^{{\square }_{\text{XY}}}\to [\tilde{\sigma }{]}_p^{{\square }_{\text{XY}}}$$

$$[\tau {]}_p^{{\square }_{\text{XZ}}}\to [\tilde{\tau }{]}_p^{{\square }_{\text{XZ}}}$$

$$\max (K_x,K_y)\ll t\ll K_z$$

$$H\to \tilde{H}$$

$$H_t\to {\tilde{H}}_t$$

$$H$$

$$H_0$$

$$H_1$$

$$H^{\text{eff}}$$

$$\mathcal{P}$$

$$\mathcal{D}$$

$$H^{\text{eff}}$$

$${\tilde{H}}_{\text{X-cube}}$$

$$-\sum _{c=\text{mancube}}{E}_{\text{mancube}}[\tilde{\sigma },\tilde{\tau }{]}_c^{\text{mancube}}$$

$$E_{\text{mancube}}$$

$$[\tilde{\sigma },\tilde{\tau }{]}_c^{\text{mancube}}$$

$$\sigma \text{-links:}$$

$$\tau \text{-links:}$$

$${\tilde{\sigma }}_{\ell }^x{\tilde{\tau }}_{\ell }^x$$

$${\tilde{\sigma }}_{\ell }^x$$

$${\tilde{\sigma }}_{\ell }^z{\tilde{\tau }}_{\ell }^z$$

$${\tilde{H}}_{\text{X-cube}}$$

$$-\sum _{\iota =+}^{\text{XY planes}}{E}_{+}[\hat{\rho }{]}_{\iota }^{{+}_{\text{XY}}}-\sum _{\iota =+}^{\text{XZ planes}}{E}_{+}[\hat{\rho }{]}_{\iota }^{{+}_{\text{XZ}}}-\sum _{c=\text{mancube}}{E}_{\text{mancube}}[\hat{\rho }{]}_c^{\text{mancube}}$$

$$\tilde{H}$$

$$\sim -E_{\square }([\tilde{\sigma }{]}_p^{{\square }_{\text{XY}}}+[\tilde{\tau }{]}_p^{{\square }_{\text{XZ}}})-E_{+}([\tilde{\sigma }{]}_{\iota }^{{+}_{\text{XY}}}+[\tilde{\tau }{]}_{\iota }^{{+}_{\text{XZ}}})$$

$$-t{\tilde{\sigma }}_{\ell }^x{\tilde{\tau }}_{\ell }^x$$

$$[\tilde{\sigma }{]}_{\iota }^{{+}_{\text{XY}}}|\Psi ⟩=[\tilde{\tau }{]}_{\iota }^{{+}_{\text{XZ}}}|\Psi ⟩=[\tilde{\sigma },\tilde{\tau }{]}_c^{\text{mancube}}|\Psi ⟩=|\Psi ⟩$$

$$[\tilde{\sigma }{]}_p^{{\square }_{\text{XY}}}$$

Full text

I Introduction

Quantum phases with fracton topological order Chamon (2005); Rasmussen et al. (2016); Xu (2006); Bravyi et al. (2011a); Haah (2011); Yoshida (2013); Vijay et al. (2015) have recently received much attention. Hsieh and Halász (2017); Prem et al. (2017); Haah (2014); Williamson (2016) These phases can be characterized as having topological excitations, called fractons Vijay et al. (2015), which are immobile (when isolated from other fractons Pretko (2017)), and for having large degeneracies (on a torus) which can be sub-extensive: e.g. $\text{degen}\sim 2^{L}$ for a 3D torus of length $L$. Fracton order has applications including self-correcting quantum memory for quantum computers Haah (2011); Bravyi and Haah (2013, 2011), and even emergent gravity given a finite density of (possibly gapless) fractons Pretko (2017).

However, most models of fracton order require many-spin interactions, which seems difficult to realize in a real material or cold atom system. For example, the so-called X-cube order requires a 12-spin interaction term Vijay et al. (2016). Recently, it has been shown Vijay (2017); Ma et al. (2017) that X-cube order can be constructed by coupling together three orthogonal stacks of intersecting two-dimensional toric codes Dennis et al. (2002), requiring only 4-spin interactions. Unfortunately, 4-spin interactions can also be difficult to realize in a physical system. On the other hand, in Ref. Kitaev (2006) Kitaev presents a honeycomb lattice model with only 2-spin interactions which, in a certain large coupling limit, has an effective Hamiltonian described by toric code. In Sec. II, we combine these ideas to produce a more realistic model of fracton order which only requires nearest-neighbor 2-spin interactions. Our model consists of two orthogonal stacks of Kitaev honeycomb models, which intersect on certain shared links. To induce fracton order, we add a large two-spin coupling on the shared links to couple together the orthogonal layers.

Ref. Vijay (2017) also showed that their isotropic coupled layer construction admits a duality describing a first order phase transition out of X-cube fracton topological order. Additionally, Ref. Vijay et al. (2016) gave dualities describing phase transitions out of an entire class of fracton phases; however, the nature of these phase transitions is currently unknown. In this work, we find two new dualities which describe two different quantum phase transitions out of the X-cube phase which both appear to be discontinuous first order transitions. The dual descriptions of the two phase transitions are described by a grid of weakly coupled 1+1D transverse-field Ising models, or a stack of decoupled 2+1D transverse-field Ising models. These dualities are summarized in Fig. 1 and are described in detail in Sec. III.

Although a single layer of a 1+1D or 2+1D transverse-field Ising model hosts a continuous quantum phase transition, this isn’t necessarily the case when they are weakly coupled together. The dual theories in our model allow the 1+1D or 2+1D Ising layers to be coupled together as long as a subdimensional Ising symmetry which flips all spins on any given layer is preserved. Therefore, a $\phi_{L}^{2}\phi_{L^{\prime}}^{2}$ term (in a stack of $\phi^{4}$ field theories) or $\sigma_{LI}^{z}\sigma_{LJ}^{z}\sigma_{L^{\prime}I}^{z}\sigma_{L^{\prime}J}^{z}$ term (in a stack of transverse-field Ising models) is symmetry allowed, where $L$ and $L^{\prime}$ denote neighboring layers, and $I$ and $J$ are nearest-neighbor sites on a layer. ($\phi_{L}\phi_{L^{\prime}}$ or $\sigma_{LI}^{z}\sigma_{L^{\prime}I}^{z}$ are odd under subdimensional symmetry and thus not allowed.) In Appendix C we give classical Monte Carlo evidence that these symmetry allowed perturbations are relevant under RG (renormalization group) and that the transitions are likely to be discontinuous first order transitions. That is, for small but finite layer coupling, we find evidence of first order transitions; but we can not rule out the possibility that smaller layer couplings could result in a continuous phase transition. We also find intermediate phases with topological (possibly with fractons) order and spontaneous translation symmetry breaking when the sign of the layer coupling is negative.

In Sec. IV we also consider an extended model where we introduce a different interaction term to couple together the intersecting honeycomb layers (and an infinitesimal perturbation to break accidental symmetry). Similar to the previous coupled layer constructions Vijay (2017); Ma et al. (2017), when this coupling is large, 3+1D $Z_{2}$ topological order is stabilized. However, when both this coupling and the coupling that resulted in X-cube order are large, a stack of 2+1D $Z_{2}$ topological order is stabilized on planes orthogonal to the original two stacks of honeycomb lattices. This result differs from the previous coupled layer constructions Vijay (2017); Ma et al. (2017) which obtain a trivial state when both couplings are large. This difference occurs because our construction is built from only two orthogonal stacks of 2+1D $Z_{2}$ topological order, while the previous constructions used three stacks. In Sec. IV.3 we give a global duality description for four (most likely first order) phase transitions in the extended model by mapping to decoupled 2+1D and 1+1D transverse-field Ising models, which is summarized in Fig. 1. In Sec. IV.4 we discuss connections to previous coupled layer constructions Vijay (2017); Ma et al. (2017). Finally, in Sec. V we conclude with possible future directions.

II Nearest-Neighbor Fracton Model

The model lives on a 3 dimensional lattice, shown in Fig. 2, of two orthogonal stacks of honeycomb lattices: one stack consisting of a honeycomb lattice on an XY plane for each $z\in\mathbb{Z}$, and the other stack on XZ planes for each $y\in\mathbb{Z}$. Some of the links (dark-blue in Fig. 2) of each honeycomb lattice will overlap with a link on an orthogonal honeycomb lattice. We will refer to these links as “shared links”, and the sites on each end of a shared link will be referred to as “shared sites”. There will be a spin-1/2 degree of freedom on every site of the honeycomb lattices, which means that there will be two spin-1/2 on each shared site. The spins on the XY and XZ planes will be labeled $\sigma^{\mu}$ and $\tau^{\mu}$, respectively.

The Hamiltonian is given by

[TABLE]

In Sec. IV, we will discuss an alternate coupling that can produce 3+1D $Z_{2}$ topological order. A rough phase diagram for this model is given in Fig. 3. The signs of all coefficients in $H$ (i.e. $K_{\mu}$ and $t$, and also $h$ and ${h^{\prime}}$ in $H^{\prime}$, which is defined later in Eq. (25)) will be taken to be positive.

$H_{K}$ is a sum of Kitaev honeycomb models Kitaev (2006) on each honeycomb lattice. $\sum_{\mu=x,y,z}$ sums over the different kinds of links ($x$=red, $y$=green, $z$=blue in Fig. 2), and $\sum_{\left<ij\right>}^{\mu\text{-links}}$ sums over all links of type $\mu$. A link is also referred to as a $\sigma$-link (or $\tau$-link) if it is a link on either a XY (or XZ) honeycomb plane. If a link is shared by both an XY and XZ honeycomb plane, then it is instead referred to as a shared-link. $K_{x}$, $K_{y}$, $K_{z}$ are couplings across the different bond types. If we only consider $H_{K}$, then the ground state will have decoupled 2+1D $Z_{2}$ topological order on each honeycomb lattice, and will be gapped when $K_{z}>K_{x}+K_{y}$ Kitaev (2006).

$H_{t}$ couples the 2D honeycomb lattices together. $H_{t}$ consists of a $\sigma^{z}_{i}\tau^{z}_{i}$ coupling on every shared site. In Sec. II.2 we argue that this term will cause loops of fluxes to condense, which will result in X-cube fracton topological order when $K_{z}$ and $t$ are large.

II.1 Decoupled Phase: XY&XZ plane 2+1D $Z_{2}$ Topological Order

Before explaining the fracton phase, we will first review the $t=0$ case with $\max(K_{x},K_{y})\ll K_{z}$. In this limit, the XY and XZ–plane hexagon lattices are completely decoupled from each other, and each form a Kitaev honeycomb model Kitaev (2006). In Fig. 4, we show a honeycomb lattice on the XY plane. When $\max(K_{x},K_{y})\ll K_{z}$, degenerate perturbation theory results in the following effective Hamiltonian (ignoring subleading terms) Kitaev (2006):

[TABLE]

[TABLE]

$\sum_{p=\square}^{\text{XY planes}}$ sums over plaquettes $p$ that lie on an XY plane and $\sum_{\iota=+}^{\text{XY plane}}$ sums over vertices $\iota$ that lie on an XY plane (see Fig. 4). $\sum_{p=\square}^{\text{XZ planes}}$ and $\sum_{\iota=+}^{\text{XZ plane}}$ are similar, but instead sum over XZ plalnes. $[\sigma]^{\square_{\text{XY}}}_{p}$ and $[\sigma]^{+_{\text{XY}}}_{\iota}$ are defined in Fig. 4 and are closely related to the plaquette/flux and vertex/charge operators in Kitaev’s toric code Dennis et al. (2002). $[\sigma]^{\square_{\text{XY}}}_{p}$ and $[\sigma]^{+_{\text{XY}}}_{\iota}$ are centered on the plaquettes $p$ and vertices $\iota$, respectively, of the XY plane rectangular lattice (gray in Fig. 4). $\tau$ operators are similar, but instead live on XZ planes. In Fig. 5, we show the 3 dimensional lattice again, but with the intersecting rectangular toric code lattices also drawn. A more detailed calculation is done in Appendix A.1, with excitation energies $E_{\square}$ and $E_{+}$ given for nonzero $t$ (Eq. (48)).

It is useful to do a change of basis and projection into the low energy Hilbert space where the energy of the large $K_{z}$ term in $H$ (Eq. (1)) is minimized with $\sigma^{z}_{i}\sigma^{z}_{j}=1$ across every $z$-link:

[TABLE]

[TABLE]

[TABLE]

Thus, we are projecting into the $\sigma^{z}_{i}\sigma^{z}_{j}=1$ subspace and doing a unitary rotation on the shared-links. $\sigma^{0}_{i}$ (e.g in $\sigma^{0}_{i}\sigma^{z}_{j}\rightarrow\tilde{\sigma}^{z}_{\ell}$) is an identity operator, and is only written to emphasize that $i$ and $j$ are the sites on the two sides of the link $\ell$. We will use a tilde to denote operators in this new basis. In this basis, the operators $[\sigma]^{\square_{\text{XY}}}_{p}\rightarrow[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ and $[\sigma]^{+_{\text{XY}}}_{\iota}\rightarrow[\tilde{\sigma}]^{+_{\text{XY}}}_{\iota}$ (also given in Fig. 4) in $\tilde{H}_{\text{toric}}$ exactly reproduce the flux and charge operators in Kitaev’s toric code Dennis et al. (2002). Therefore, when $t=0$ and $\max(K_{x},K_{y})\ll K_{z}$, our model (Eq. (3)) is in a decoupled phase where every honeycomb lattice on an XY or XZ plane hosts its own gapped 2+1D $Z_{2}$ topological order. This is the “XY&XZ plane 2+1D $Z_{2}$ QSL” phase in Fig. 3.

II.2 X-Cube Fracton Phase

In this section, we will explain how X-cube fracton topological order results from our model $H$ (Eq. (1)) in the limit

[TABLE]

More generally, in the Appendix A.2 we show that $\max(K_{x},K_{y})\ll\min(t,K_{z})$ is a sufficient condition for X-cube order. However, Eq. (7) allows us to apply the results from the previous section and is more closely connected to the previous layer constructions in Ref. Vijay (2017); Ma et al. (2017), and so we will discuss this limit in the main text and consider the more general condition (Eq. (50)) in the appendix.

Thus, we will start from the decoupled $t=0$ limit discussed in the previous section, and then make $t$ very large. In the toric code basis (Eq. (5)), $H$ (Eq. (1)) becomes:

[TABLE]

where $\tilde{H}_{\text{toric}}$ was given in Eq. (6). We see that $\tilde{H}_{t}$ couples the XY and XZ–plane toric codes together on the shared-links (gray links parallel to x-axis in Fig. 5) with a $\tilde{\sigma}^{x}_{\ell}\tilde{\tau}^{x}_{\ell}$ coupling. When acting on the decoupled $t=0$ ground state, $\tilde{\sigma}^{x}_{\ell}\tilde{\tau}^{x}_{\ell}$ creates two $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ and two $[\tilde{\tau}]^{\square_{\text{XZ}}}_{p}$ excitations on the four plaquettes (of the gray cubic lattice in Fig. 5) that neighbor the link $\ell$. When $t$ is large, these excitations condense and the result is X-cube fracton topological order. (A similar condensation process is described in Ref. Vijay (2017) and Ma et al. (2017); see Sec. IV.4 for connections.)

We will derive this result using degenerate perturbation theory (DPT). First, we split $H$ (Eq. (1)) into

[TABLE]

where $H_{0}$ is the unperturbed Hamiltonian, $H_{1}$ is the perturbing Hamiltonian, and $\tilde{H}_{\text{toric}}$ (Eq. (6)) was derived in the previous section. Eq. (7) is sufficient to ensure that DPT converges.

The effective Hamiltonian of the low energy Hilbert space is given by (see Appendix B for details):

[TABLE]

$\mathcal{P}$ projects into the low energy Hilbert space where $H_{0}=\tilde{H}_{t}$ is minimized. $\mathcal{D}$ projects onto states not in the low energy Hilbert space and penalizes energy excitations by a factor $1/(H_{0}-E_{0})$ where $E_{0}$ is the (degenerate) ground state energy of $H_{0}$. “$\cdots$” contains many terms (e.g. $-\frac{1}{2}\mathcal{P}H_{1}\mathcal{D}^{2}H_{1}\mathcal{P}H_{1}\mathcal{P}$); but these more complicated terms will be subleading corrections in this section. $[\tilde{\sigma}]^{+_{\text{XY}}}_{\iota}$ and $[\tilde{\tau}]^{+_{\text{XZ}}}_{\iota}$ (in $H_{1}=\tilde{H}_{\text{toric}}$) commute with $H_{t}$ (Eq. (3)), which implies that they do not create excitations out of the low energy Hilbert space. However, $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ and $[\tilde{\tau}]^{\square_{\text{XZ}}}_{p}$ do not commute with $\tilde{H}_{t}$, since $\tilde{H}_{t}$ creates two $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ and two $[\tilde{\tau}]^{\square_{\text{XZ}}}_{p}$ excitations. Viewed differently, $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ and $[\tilde{\tau}]^{\square_{\text{XZ}}}_{p}$ create $\tilde{H}_{t}$ excitations, which are not in the low energy Hilbert space. Therefore, $[\tilde{\sigma}]^{+_{\text{XY}}}_{\iota}$ and $[\tilde{\tau}]^{+_{\text{XZ}}}_{\iota}$ will be contained in $H^{\text{eff}}$ but $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ and $[\tilde{\tau}]^{\square_{\text{XZ}}}_{p}$ will not.

We must now consider higher order terms ($n\geq 2$) in $H^{\text{eff}}$ (Eq. (11)) to discover what role $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ and $[\tilde{\tau}]^{\square_{\text{XZ}}}_{p}$ end up playing. The first nonzero term involving $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ or $[\tilde{\tau}]^{\square_{\text{XZ}}}_{p}$ is $E_{\text{\mancube}}[\tilde{\sigma},\tilde{\tau}]^{\text{\mancube}}_{c}$, which is obtained at fourth order ($n=4$), and the resulting effective Hamiltonian is

[TABLE]

$\sum_{c=\text{\mancube}}$ sums over the cubes $c$ shown in Fig. 5. $[\tilde{\sigma},\tilde{\tau}]^{\text{\mancube}}_{c}$ is the product of two XY-plane plaquette $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ and two XZ-plane plaquette $[\tilde{\tau}]^{\square_{\text{XZ}}}_{p}$ operators around a cube (Fig. 6). Note that $[\tilde{\sigma},\tilde{\tau}]^{\text{\mancube}}_{c}$ commutes with $\tilde{H}_{t}$, and can therefore appear in the effective Hamiltonian since it preserves the low energy Hilbert space. The “$\cdots$” in $H^{\text{eff}}$ denotes other subleading terms generated by DPT which are not important. Another round of DPT on $H^{\text{eff}}$ with an unperturbed Hamiltonian $H_{0}^{\prime}=\tilde{H}_{\text{X-cube}}$ (and perturbation $H_{1}^{\prime}=H^{\text{eff}}-\tilde{H}_{\text{X-cube}}$) will eliminate terms in “$\cdots$” which anticommute with $\tilde{H}_{\text{X-cube}}$, resulting in an exactly solvable effective Hamiltonian with the same ground state as $\tilde{H}_{\text{X-cube}}$. The “$\sim$” in the energy expressions means that these expressions are only asymptotically correct (in the $\max(K_{x},K_{y})\ll t\ll K_{z}$ limit); e.g. factors of two are ignored.

It is again useful to do a change of basis and projection into the low energy Hilbert space where now the energy of $H_{0}$ is minimized with $\tilde{\sigma}^{x}_{\ell}\tilde{\tau}^{x}_{\ell}=1$ on every shared-link:

[TABLE]

[TABLE]

[TABLE]

Thus, we are relabeling the $\sigma$ and $\tau$–links, but on the shared-links we are projecting into the $\tilde{\sigma}^{x}_{\ell}\tilde{\tau}^{x}_{\ell}=1$ subspace. A hat over operators is used to denote operators in this new basis. See Fig. 6 and Fig. 7 for descriptions of the operators in $\hat{H}_{\text{X-cube}}$. In this basis, $\hat{H}_{\text{X-cube}}$ almost exactly reproduces the X-cube model introduced in Vijay et al. (2016). The only difference is that $\hat{H}_{\text{X-cube}}$ does not include the $[\hat{\rho}]^{+_{\text{YZ}}}_{\iota}\equiv[\hat{\rho}]^{+_{\text{XY}}}_{\iota}[\hat{\rho}]^{+_{\text{XZ}}}_{\iota}$ operator. However, this does not change the phase of the Hamiltonian; it merely adds anisotropy. Only two of $[\hat{\rho}]^{+_{\text{XY}}}_{\iota}$, $[\hat{\rho}]^{+_{\text{XZ}}}_{\iota}$, and $[\hat{\rho}]^{+_{\text{YZ}}}_{\iota}$ are necessary for X-cube order, since the third is a product of the other two. Therefore, when $\max(K_{x},K_{y})\ll\min(K_{z},t)$, our model (Eq. (1)) has X-cube fracton topological order, as claimed in the phase diagram Fig. 3. Note that $\hat{H}_{\text{X-cube}}$ is actually a rather simple Hamiltonian in the sense that all of its terms commute. Its ground state is specified by $[\hat{\rho}]^{+_{\text{XY}}}_{\iota}=[\hat{\rho}]^{+_{\text{XZ}}}_{\iota}=[\hat{\rho}]^{\text{\mancube}}_{c}=1$ and has degeneracy $2^{6L-3}$ on a 3D torus Vijay et al. (2016); Ma et al. (2017). It is also interesting to note that the X-cube operators generated by DPT ($[\sigma,\tau]^{\text{\mancube}}_{p}$, $[\sigma]^{+_{\text{XY}}}_{\iota}$ and $[\tau]^{+_{\text{XZ}}}_{\iota}$) all commute with the original Hamiltonian (Eq. (1)).

$[\hat{\rho}]^{\text{\mancube}}_{c}$ is the fracton number operator. Its excitations are of particular interest because they can only be created in groups of four Vijay et al. (2015, 2016). This can be seen from Fig. 6 where the application of a single $\hat{\rho}^{x}_{\ell}$ operator on any link $\ell$ will excite four $[\hat{\rho}]^{\text{\mancube}}_{c}$ fracton operators on the four cubes that share the link $\ell$. If we have an isolated fracton, then due to the $Z_{2}$ nature of the fracton operator, the application of a $\hat{\rho}^{x}_{\ell}$ operator will annihilate the isolated fracton but create three other fractons. Therefore, an isolated fracton can not move without creating additional excitations. An isolated fracton is immobile; only a pair of fractons can move without creating additional excitations.

III Duality and Phase Transitions out of X-cube Order

III.1 2+1D Duality

In this subsection we will discuss a duality (similar to the one in Ref. Vijay (2017), see Sec. IV.4 for connections) which describes a (likely first order) phase transition from X-cube fracton order to decoupled XY&XZ plane 2+1D $Z_{2}$ topological order. The duality will be applied to $\tilde{H}$ (Eq. (8)):

[TABLE]

The last line above is just a short-hand representation of $\tilde{H}$, and is only shown for convenience.

Recall that $\tilde{H}_{\text{toric}}$ (Eq. (6)) describes decoupled stacks of 2+1D toric codes along the XY and XZ planes. (The tilde over operators is still used to distinguish from the Pauli operators in the original Hamiltonian (Eq. (1)).) The XY-plane toric codes are composed of Pauli operators $\tilde{\sigma}^{\mu}_{\ell}$ on the links of the XY-plane square lattices, and are described by plaquette/flux operators $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ and vertex/charge operators $[\tilde{\sigma}]^{+_{\text{XY}}}_{\iota}$ (Fig. 4). The XZ-plane toric codes are similar, but are defined using $\tilde{\tau}$ instead of $\tilde{\sigma}$. The XY and XZ–plane toric code square lattices intersect on the links parallel to the x-axis, which we have previously referred to as shared-links. $\tilde{H}_{t}$ induces X-cube order at large $t$ by coupling these shared-links together via a $t\,\tilde{\sigma}^{x}_{\ell}\tilde{\tau}^{x}_{\ell}$ coupling.

First, note that that $[\tilde{\sigma}]^{+_{\text{XY}}}_{\iota}$, $[\tilde{\tau}]^{+_{\text{XZ}}}_{\iota}$, and $[\tilde{\sigma},\tilde{\tau}]^{\text{\mancube}}_{c}$ all commute with $\tilde{H}$. Recall that $[\tilde{\sigma},\tilde{\tau}]^{\text{\mancube}}_{c}$ is the product of two XY-plane and two XZ-plane flux operators ($[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ and $[\tilde{\tau}]^{\square_{\text{XZ}}}_{p}$) that boarder the cube $c$. (Note, there are no YZ-plane flux operators on the other two faces of the cube.) Therefore, we can work in the restricted Hilbert space where

[TABLE]

This is similar to imposing a charge-free gauge constraint when one considers dualities of pure $Z_{2}$ gauge theory. For example, our restricted Hilbert space does not allow any fracton excitations, which are excitations of $[\tilde{\sigma},\tilde{\tau}]^{\text{\mancube}}_{c}$.

The dual theory will consist of a stack of decoupled YZ-plane square lattices with vertices $\hat{\iota}$ at the center of the x-axis links $\ell$ of the original cubic lattice. We introduce a new Pauli operator $\eta^{\mu}_{\hat{\iota}}$ which lives on the dual lattice vertices (Fig. 1 (b)). $\eta^{\mu}_{\hat{\iota}}$ will only be coupled along YZ planes; different YZ planes will be completely decoupled from each other. The duality is given by

[TABLE]

Therefore, the flux (XY plane $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ or XZ plane $[\tilde{\tau}]^{\square_{\text{XZ}}}_{p}$) through the plaquette $p$ is mapped to $\eta^{z}_{\hat{\iota}}\eta^{z}_{\hat{\iota}^{\prime}}$ where $\hat{\iota}$ and $\hat{\iota}^{\prime}$ reside at the center of the two x-axis links of the plaquette $p$. The coupling $\tilde{\sigma}^{x}_{\ell}\tilde{\tau}^{x}_{\ell}$ (on an x-axis link) is mapped to $\eta^{x}_{\hat{\iota}}$, which resides at the center of the link $\ell$. $[\tilde{\sigma}]^{+_{\text{XY}}}_{\iota}$, $[\tilde{\tau}]^{+_{\text{XZ}}}_{\iota}$, and $[\tilde{\sigma},\tilde{\tau}]^{\text{\mancube}}_{c}$ are mapped to the identity, which solves the Hilbert space constraint Eq. (17).

Note that the algebra (e.g. commutation relations) of the operators on the left hand side (LHS) and right hand side (RHS) of Eq. (18) are the same. Furthermore, the LHS and RHS both have exactly one effective qubit per unit cell. The LHS originally has 4 qubits (two $\tilde{\sigma}$ and two $\tilde{\tau}$), but this gets subtracted by 3 qubits due to the 3 restricted Hilbert space conditions (Eq. (17)). Therefore, the transformation (Eq. (18)) is an algebra isomorphism (on an infinite lattice) between the algebra of $\tilde{\sigma}^{\mu}_{\ell}$ and $\tilde{\tau}^{\mu}_{\ell}$ modulo the Hilbert space restriction (Eq. (17)), and the algebra of $\eta^{\mu}_{\hat{\iota}}$. Thus, Eq. (18) is a valid duality transformation.

The dual Hamiltonian is described by decoupled 2+1D transverse field Ising models on the dual YZ-plane square lattices:

[TABLE]

$\sum_{\langle\hat{\iota}\hat{\iota}^{\prime}\rangle}^{\text{YZ planes}}$ sums over nearest-neighbor dual lattice vertices that lie on the same YZ plane. We see that the decoupled phase (XY & XZ plane $Z_{2}$ topological order) with small coupling $t$ is dual to the ordered phase of $\tilde{H}^{\text{dual}}$, and the X-cube phase with large $t$ is dual to the disordered phase of $\tilde{H}^{\text{dual}}$.

In the original model $\tilde{H}$ (Eq. (16)), the $t$ coupling, which couples together orthogonal layers of toric code, creates a YZ-plane loop of four fluxes (two XY and two XZ–plane fluxes: Fig. 1 (b)). In the previous subsection (Sec. II.2) we use degenerate perturbation theory to show that at strong coupling, this YZ-plane composite flux loop condenses, and X-cube order results. Similar to the explanations in Ref. Vijay (2017); Ma et al. (2017), this occurs because the XY and XZ–plane toric code charge excitations ($[\tilde{\sigma}]^{+_{\text{XY}}}_{\iota}$ and $[\tilde{\tau}]^{+_{\text{XZ}}}_{\iota}$) can no longer move in the x direction since they have nontrivial braiding statistics with the YZ-plane composite flux loop. Only a bound state of an XY and an XZ–plane toric code charge excitation can move in the x direction; but this composite particle can no longer move in any other direction. One can also consider a pair of open YZ-plane flux loop excitations, which have immobile fracton excitations $[\tilde{\sigma},\tilde{\tau}]^{\text{\mancube}}_{c}$ (Eq. (6)) at the open endpoints. Before the composite flux condensation, open flux loops have energy proportional to their length; however, after condensation, the bulk of the open flux loop is indistinguishable from the bulk of a closed flux loop, which carries no energy. Thus, arbitrarily long open YZ-plane flux loops only carry a finite energy cost at their boundaries. We see that in $\tilde{H}^{\text{dual}}$, the $t$ coupling flips a $\eta^{z}$ spin, which will either create or deform a domain wall on a YZ plane. Thus, the condensation of YZ-plane flux loops is dual to a proliferation of YZ-plane domain walls in 2+1D YZ-plane Ising models.

Although the spontaneous symmetry breaking phase transition of a single layer of the 2+1D transverse-field Ising model is continuous, the same transition for a stack of Ising models (i.e. $\tilde{H}^{\text{dual}}$) appears to be first order once they are weakly coupled. This is because a $\phi_{L}^{2}\phi_{L^{\prime}}^{2}$ term (in a stack of $\phi^{4}$ field theories) or $\sigma_{LI}^{z}\sigma_{LJ}^{z}\sigma_{L^{\prime}I}^{z}\sigma_{L^{\prime}J}^{z}$ term (in a stack of transverse-field Ising models) is allowed by the subsystem symmetry, where $L$ and $L^{\prime}$ denote different layers, and $I$ and $J$ are nearest-neighbor sites on a layer. ($\phi_{L}\phi_{L^{\prime}}$ or $\sigma_{LI}^{z}\sigma_{L^{\prime}I}^{z}$ are odd under subdimensional symmetry and thus not allowed.) $\phi^{2}$ has scaling dimension $\sim 1.4$, and thus $\phi_{L}^{2}\phi_{L^{\prime}}^{2}$ can be expected to have dimension $\sim 1.4\times 2\sim 2.8$, which is less than the dimension (2+1) of the critical Ising layers. Thus $\phi_{L}^{2}\phi_{L^{\prime}}^{2}$ is likely to be relevant under RG at the critical point, which implies that the transition of weakly coupled 2+1D Ising layers is not described by the decoupled transition. According to the duality mapping (Eq. (18)), the relevant operator corresponds to a product of two $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ or $[\tilde{\tau}]^{\square_{\text{XZ}}}_{p}$ operators. In Appendix C we provide classical Monte Carlo evidence that these symmetry allowed perturbations are relevant under RG, and that this transition and its dual are likely to be discontinuous and first order.

III.2 1+1D Duality

In this subsection we will discuss a duality describing a (likely first order) phase transition from the X-cube phase to a phase of decoupled YZ-plane 2+1D $Z_{2}$ topological order. To do this, we will consider the X-cube model, but with a dimension-1 particle hopping term (and a perturbation):

[TABLE]

where $\sum_{\ell}^{\text{x-axis links}}$ sums over all links $\ell$ which are parallel to the x-axis. The hopping term $h$ (which we also consider in Sec. IV) hops the $[\hat{\rho}]^{+_{\text{XY}}}_{\iota}$ and $[\hat{\rho}]^{+_{\text{XZ}}}_{\iota}$ excitations (Fig. 7) along the x-axis, which will condense the X-cube ordered phase into YZ-plane 2+1D $Z_{2}$ topological order. A small $E^{\text{YZ}}_{\square}\,[\hat{\rho}]^{\square_{\text{YZ}}}_{p}$ perturbation is included to prevent fine tuning, where $[\hat{\rho}]^{\square_{\text{YZ}}}_{p}$ is a YZ-plane toric code flux operator (Fig. 9). Similar to $H_{h^{\prime}}$ (Eq. (27)) in a later section, without the $E^{\text{YZ}}_{\square}$ term, the large $h$ phase would have an unphysical degeneracy $\sim 2^{L^{2}}$ due to accidental symmetries which prevent a $[\hat{\rho}]^{\square_{\text{YZ}}}_{p}$ flux term from being generated. Terms other than $E^{\text{YZ}}_{\square}\,[\hat{\rho}]^{\square_{\text{YZ}}}_{p}$ can also work; however, this is the simplest choice for understanding the duality.

$\hat{H}_{\text{X-cube}}$ describes the X-cube model Vijay et al. (2016) with $\hat{\rho}^{\mu}$ operators on the links of a cubic lattice. (The hat over operators is used to distinguish them from the Pauli operators in the original Hamiltonian (Eq. (1) and previous section).) $[\hat{\rho}]^{\text{\mancube}}_{c}$ (Fig. 6) measures the fracton number and is a product of the twelve $\hat{\rho}^{z}_{\ell}$ operators on the edges of the cube $c$. $[\hat{\rho}]^{+_{\text{XY}}}_{\iota}$, $[\hat{\rho}]^{+_{\text{XZ}}}_{\iota}$, and $[\hat{\rho}]^{+_{\text{YZ}}}_{\iota}=[\hat{\rho}]^{+_{\text{XY}}}_{\iota}[\hat{\rho}]^{+_{\text{XZ}}}_{\iota}$ (Fig. 7) are products of four $\hat{\rho}^{x}_{\ell}$ operators and correspond to dimension-1 particle excitations, which can only move along 1-dimensional lines. (The absence of $[\hat{\rho}]^{+_{\text{YZ}}}_{\iota}$ in $\hat{H}_{\text{X-cube}}$ is not important.)

$[\hat{\rho}]^{+_{\text{YZ}}}_{\iota}$ and $[\hat{\rho}]^{\text{\mancube}}_{c}$ commute with $\hat{H}_{\text{hop}}$ (Eq. (20)), and therefore we can work in a restricted Hilbert space where

[TABLE]

The dual theory will consist of a grid of weakly coupled Ising chains along the x-axis with the same lattice sites $\iota$ as the original cubic lattice. We introduce a new Pauli operator $\mu^{\alpha}_{\iota}$ which lives on the vertices of the cubic lattice. The duality is given by

[TABLE]

Thus, the dimension-1 particle ($[\hat{\rho}]^{+_{\text{XY}}}_{\iota}$ and $[\hat{\rho}]^{+_{\text{XZ}}}_{\iota}$) number operator is mapped to a spin flip operator $\mu^{x}_{\iota}$; and the dimension-1 particle hopping $\hat{\rho}^{z}_{\ell}$ is mapped to a domain wall number operator $\mu^{z}_{\iota}\mu^{z}_{\iota^{\prime}}$. $[\hat{\rho}]^{+_{\text{YZ}}}_{\iota}$ and $[\tilde{\sigma},\tilde{\tau}]^{\text{\mancube}}_{c}$ are mapped to the identity, which solves the Hilbert space constraint Eq. (22). The YZ-plane flux operator $[\hat{\rho}]^{\square_{\text{YZ}}}_{p}$ is mapped to a product of four $\mu^{z}_{\iota}$ operators at the corners of the plaquette $p$. The duality mapping forms an algebra isomorphism (on an infinite lattice) of the left hand side (LHS) of Eq. (23) (modulo the Hilbert space restriction (Eq. (22))), and the right hand side (RHS). For example, both the LHS and RHS have one effective qubit per unit cell. The LHS originally has 3 qubits ($\hat{\rho}^{\mu}_{\ell}$ along the three different directions), but this gets subtracted by 2 qubits due to the 2 Hilbert space conditions (Eq. (22)).

The dual Hamiltonian is described by a grid of 1+1D transverse field Ising models parallel to the x-axis (Fig. 1 (c)):

[TABLE]

We see that the X-cube phase with small hopping $h$ is dual to the disordered phase of $\hat{H}_{\text{hop}}^{\text{dual}}$, and the YZ-plane topological order phase with large $h$ is dual to the ordered phase of $\hat{H}_{\text{hop}}^{\text{dual}}$.

In the original model $\hat{H}_{\text{hop}}$ (Eq. (20)), the $h$ hopping hops dimension-1 particles (a composite excitation of $[\hat{\rho}]^{+_{\text{XY}}}_{\iota}$ and $[\hat{\rho}]^{+_{\text{XZ}}}_{\iota}$) parallel to the x-axis. For large hopping strength $h$, these x-axis particles condense, and YZ-plane 2+1D $Z_{2}$ topological order results. This is shown in a latter section (Sec. IV.2) using degenerate perturbation theory with a more general model including a $t$ term, which is taken to be zero in this section. We now give a physical explanation. The x-axis particle condensation confines fracton excitations which are created by YZ-plane membrane operators, since these membrane operators have nontrivial statistics with the x-axis particle. This condensation also allows the y-axis particle $[\hat{\rho}]^{+_{\text{XY}}}_{\iota}$ to move along the z direction, since it can do this by splitting up into a z-axis particle (which can move in the z direction) and an x-axis particle, which costs no energy after its condensation. Thus, the y-axis and z-axis particles become YZ-plane toric code charge excitations. Therefore, after the x-axis particle condensation, there are no dimension-1 particles, which implies that the fracton excitations created by XY and XZ–plane membrane operators are no longer finite energy excitations. To understand this, note that in the X-cube phase, there is no excitation (with a nontrivial braiding statistic) that can be braided through the membrane along a small closed loop. However, after the x-axis particles are condensed and YZ-plane toric code charge excitations are permitted, the YZ-plane charge excitations can braid through an XY-plane membrane. As a result, an XY-plane (or XZ-plane) membrane will have an energy cost proportional to the number of YZ-plane charges that can braid around it, which is proportional to the length of the membrane in the x direction. Mathematically, this energy cost results from the YZ-plane toric code flux term $E^{\text{YZ}}_{\square}\,[\hat{\rho}]^{\square_{\text{YZ}}}_{p}$ in $\hat{H}_{\text{hop}}$ (Eq. (20)). (The same physics explains why flux loops in 3+1D toric code have energy cost proportional to their length.) However, an XY-plane (or XZ-plane) membrane operator with small length (e.g. =1) in the x direction will have a finite energy cost; these operators create composite fracton excitations which are free to move along a YZ plane. After the x-axis particles condense, this composite fracton excitation becomes the YZ-plane toric code flux excitation.

In the dual theory $\hat{H}_{\text{hop}}^{\text{dual}}$, $h$ hops excitations of the spin flip operator $\mu^{x}_{\iota}$, which spontaneously breaks $Z_{2}$ Ising symmetry when condensed since these excitations are conserved modulo two due to the $Z_{2}$ Ising symmetry. The $E^{\text{YZ}}_{\square}$ term couples four neighboring 1+1D chains together, which is a marginal under RG operator at the Ising phase transition. This can be understood by mapping the Ising chains to free Majorana chains using a Jordan-Wigner transformation; the $E^{\text{YZ}}_{\square}$ term is then mapped to a mass-mass term $(\chi_{l}^{T}\sigma^{2}\chi_{l})(\chi_{l^{\prime}}^{T}\sigma^{2}\chi_{l^{\prime}})$. Since $\chi_{l}^{T}\sigma^{2}\chi_{l}$ has scaling dimension $1$, the $E^{\text{YZ}}_{\square}$ term has scaling dimension $2\times 1$, which is equal to the number of space-time dimensions of each chain (1+1). Thus, the $E^{\text{YZ}}_{\square}$ term is marginal under RG at the fixed point of critical decoupled Ising chains. In Appendix C we show classical Monte Carlo evidence that these symmetry allowed perturbations are likely to be relevant under RG, and that this transition and its dual are likely to be discontinuous and first order.

In the next section, we introduce an extended honeycomb-based model $H^{\prime}$ (Eq. (25)) which also exhibits this phase transition with a similar duality description at low energies.

IV Extended Model with 3+1D $Z_{2}$ Topological Order

To induce 3+1D $Z_{2}$ topological order instead of fracton order, we can consider adding to $H=H_{K}+H_{t}$ (Eq. (1)) a coupling $H_{h}$ and a small perturbation, for which we will choose $H_{{h^{\prime}}}$ as an example:

[TABLE]

A rough phase diagram for this model is given in Fig. 8. See Sec. IV.4 for connections to previous work.

$H_{h}$ is an alternative coupling of the 2D honeycomb lattices. $\sum_{\left<ij\right>}^{\text{shared-links}}$ sums over all shared links. As we will explain in Sec. IV.1, this term will cause a composite charge excitation to condense, resulting in a phase which is fine tuned (with unphysical degeneracy $\sim 2^{L^{2}}$ on a torus). When generic perturbations are added (e.g. ${h^{\prime}}$), the resulting phase is 3+1D $Z_{2}$ topological order.

$H_{{h^{\prime}}}$ will be considered as an example of a generic perturbation to break the unphysical degeneracy mentioned in the previous paragraph. The particular choice of $H_{{h^{\prime}}}$ is not important.

In the tilde basis (which includes a low energy Hilbert space projection from degenerate perturbation theory), $H^{\prime}$ is given by (ignoring subleading terms)

[TABLE]

$\tilde{H}_{\text{toric}}$ (Eq. (6)) and $\tilde{H}_{t}$ (Eq. (9)) were discussed previously. We see that $\tilde{H}_{h}$ couples the XY and XZ–plane toric codes together on the shared-links with a $\tilde{\sigma}^{z}_{\ell}\tilde{\tau}^{z}_{\ell}$ coupling. If $\ell=\langle\iota\iota^{\prime}\rangle$ is a cubic lattice link between vertices $\iota$ and $\iota^{\prime}$, then when acting on the decoupled $t=h=0$ ground state, $\tilde{\sigma}^{z}_{\ell}\tilde{\tau}^{z}_{\ell}$ creates the following four excitations: $[\tilde{\sigma}]^{+_{\text{XY}}}_{\iota}$, $[\tilde{\sigma}]^{+_{\text{XY}}}_{\iota^{\prime}}$, $[\tilde{\tau}]^{+_{\text{XZ}}}_{\iota}$, $[\tilde{\tau}]^{+_{\text{XZ}}}_{\iota^{\prime}}$. When $t$ is small and $h$ is large (and ${h^{\prime}}$ is small but finite), these excitations condense and the result is 3+1D $Z_{2}$ topological order.

IV.1 3+1D $Z_{2}$ Topological Order

In this section, we will explain how 3+1D $Z_{2}$ topological order results from $H^{\prime}$ (Eq. (25)) in the limit

[TABLE]

and with infinitesimal ${h^{\prime}}$ (to break accidental degeneracies). In Appendix A.3 we will consider more general limits. However, $\max(K_{x},K_{y})\ll h\ll K_{z}$ allows us to apply the results from Sec. II.1 and is closely connected to the previous layer constructions in Ref. Vijay (2017); Ma et al. (2017), and so we will discuss this limit here.

Again, we will derive this result using degenerate perturbation theory (DPT). First, we split $H^{\prime}$ (Eq. (25)) into

[TABLE]

where $\tilde{H}_{\text{toric}}$ (Eq. (6)) describes the decoupled toric codes on XY and XZ planes. Eq. (31) is sufficient to ensure that DPT converges. $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ and $[\tilde{\tau}]^{\square_{\text{XZ}}}_{p}$ (in $\tilde{H}_{\text{toric}}$) commute with $H_{0}=\tilde{H}_{h}$, and will therefore appear in the effective Hamiltonian $H^{\text{eff}}$. These operators will play the role of flux operators on XY and XZ planes in 3+1D toric code. However, $[\tilde{\sigma}]^{+_{\text{XY}}}_{\iota}$ and $[\tilde{\tau}]^{+_{\text{XZ}}}_{\iota}$ will not appear in $H^{\text{eff}}$ since they do not commute with $H_{0}$, which means that they create excitations out of the low energy Hilbert space. However, their product $[\tilde{\sigma},\tilde{\tau}]^{*}_{\iota}\equiv[\tilde{\sigma}]^{+_{\text{XY}}}_{\iota}[\tilde{\tau}]^{+_{\text{XZ}}}_{\iota}$ does commute with $H_{0}$ and will indeed appear in $H^{\text{eff}}$. In a new basis, $[\tilde{\sigma},\tilde{\tau}]^{*}_{\iota}$ will take the form of a charge operator in 3+1D toric code. Thus, DPT will produce the following effective Hamiltonian

[TABLE]

(Actually, two applications of DPT are necessary to generate $E^{\text{YZ}}_{\square}\,[\tilde{\sigma},\tilde{\tau}]^{\square_{\text{YZ}}}_{p}$, as is done in Appendix A.3; this term will be explained soon.) Again, we perform a change of basis and projection into the low energy Hilbert space where now the energy of $H_{0}$ is minimized with $\tilde{\sigma}^{z}_{\ell}\tilde{\tau}^{z}_{\ell}=1$ on every shared-link:

[TABLE]

[TABLE]

A hat over operators is also used to denote operators in this basis. In this basis, $\hat{H}_{\text{3D toric}}$ has the exact same form as 3+1D toric code, which has 3+1D $Z_{2}$ topological order.

However, we have not yet discussed the $E^{\text{YZ}}_{\square}[\tilde{\sigma},\tilde{\tau}]^{\square_{\text{YZ}}}_{p}$ term in $\tilde{H}_{\text{3D toric}}$. In the new basis (Eq. (35)), $[\hat{\rho}]^{\square_{\text{YZ}}}_{p}$ takes the form of a flux operator on YZ planes in 3+1D toric code. This term is necessary in 3+1D toric code; without it, $\hat{H}_{\text{3D toric}}$ would have an unphysical degeneracy $\sim 2^{L^{2}}$ on a 3D torus. But without $H_{{h^{\prime}}}$, $[\tilde{\sigma},\tilde{\tau}]^{\square_{\text{YZ}}}_{p}$ can not be generated by DPT (due to accidental symmetries). Therefore, we add an infinitesimal perturbation $H_{{h^{\prime}}}$ to break this degeneracy. The choice of perturbation is not important, and we only choose $H_{{h^{\prime}}}$ as a simple example. Most sufficiently general and sufficiently small perturbations will generate $[\tilde{\sigma},\tilde{\tau}]^{\square_{\text{YZ}}}_{p}$ and not affect the resulting phase.

IV.2 YZ Plane 2+1D $Z_{2}$ Topological Order

In this section, we will explain why 2+1D $Z_{2}$ topological order on YZ planes results from $H^{\prime}$ (Eq. (25)) in the limit

[TABLE]

with infinitesimal ${h^{\prime}}$ (to break accidental degeneracies). This phase is the same as a stack of decoupled toric codes on YZ planes. In the Appendix A.4 we will consider limits more general than Eq. (37). Starting again from the decoupled $t=h=0$ limit discussed in Sec. II.1, we make both $h$ and $t$ very large and perform DPT on $\tilde{H}^{\prime}$ (Eq. (28)). $\tilde{H}_{t}$ (Eq. (9)) and $\tilde{H}_{h}$ (Eq. (29)) couple the XY and XZ–plane toric codes together on the shared-links with $\tilde{\sigma}^{x}_{\ell}\tilde{\tau}^{x}_{\ell}$ and $\tilde{\sigma}^{z}_{\ell}\tilde{\tau}^{z}_{\ell}$ couplings. When acting on the decoupled $t=h=0$ ground state, these couplings create toric code charge and flux excitations, and when $t$ and $h$ are large, these excitations condense and the result is 2+1D $Z_{2}$ topological order on YZ planes.

To understand this, we again apply degenerate perturbation theory (DPT) by splitting $H^{\prime}$ (Eq. (25)) into

[TABLE]

Eq. (37) is sufficient to ensure that DPT converges. No operator in $\tilde{H}_{\text{toric}}$ (Eq. (6)) will appear in $H^{\text{eff}}$ since they do not commute with $H_{0}$; i.e. they create excitations out of the low energy Hilbert space. However, $[\tilde{\sigma},\tilde{\tau}]^{\text{\mancube}}_{c}$ and $[\tilde{\sigma},\tilde{\tau}]^{+_{\text{YZ}}}_{\iota}$ ($=[\tilde{\sigma},\tilde{\tau}]^{*}_{\iota}$ from Eq. (34)) are generated by DPT, and the effective Hamiltonian is

[TABLE]

(Actually, two applications of DPT are necessary to generate $E^{\text{YZ}}_{+}\,[\tilde{\sigma},\tilde{\tau}]^{\square_{\text{YZ}}}_{p}$, as is done in Appendix A.4.) We perform a change of variables and projection into the low energy Hilbert space with $\tilde{\sigma}^{x}_{\ell}\tilde{\tau}^{x}_{\ell}=\tilde{\sigma}^{z}_{\ell}\tilde{\tau}^{z}_{\ell}=1$ on every shared-link:

[TABLE]

[TABLE]

A hat over operators is also used to denote operators in this basis. In this basis, $\hat{H}^{\text{YZ}}_{\text{toric}}$ has the exact same form as decoupled 2+1D toric codes on YZ planes. $[\hat{\rho}]^{+_{\text{YZ}}}_{\iota}$ and $[\hat{\rho}]^{\square_{\text{YZ}}}_{p}$ are charge and flux operators on YZ planes in 2+1D toric code, respectively.

But just like $[\tilde{\sigma},\tilde{\tau}]^{\square_{\text{YZ}}}_{p}$ (Eq. (33)) in the previous section, $[\tilde{\sigma},\tilde{\tau}]^{\square_{\text{YZ}}}_{p}$ in Eq. (39) can not be generated from $H^{\prime}$ (Eq. (25)) with ${h^{\prime}}=0$. However, $H^{\prime}$ with ${h^{\prime}}=0$ does generate $[\tilde{\sigma},\tilde{\tau}]^{\text{\mancube}}_{c}$, which corresponds to a product of flux operators on neighboring YZ planes in the new basis (Eq. (42)); that is, $[\tilde{\sigma},\tilde{\tau}]^{\text{\mancube}}_{c}\rightarrow[\hat{\rho}]^{\text{\mancube}}_{c}=[\hat{\rho}]^{\square_{\text{YZ}}}_{p}[\hat{\rho}]^{\square_{\text{YZ}}}_{p^{\prime}}$ where YZ-plane plaquettes $p$ and $p^{\prime}$ lie on the boundary of the cube $c$. Without perturbations such as $H_{h^{\prime}}$, $H^{\text{eff}}$ in Eq. (39) (and therefore also $H^{\prime}$) would have an unphysical degeneracy $\sim 2^{L^{2}}$ on a torus. Again, the choice of perturbation $H_{{h^{\prime}}}$ is not important.

IV.3 Global Duality of the Extended Model

In this subsection we discuss a global duality which describes the four (likely first order) phase transitions along the edge of Fig. 8. To do this, we merge $\tilde{H}$ (Eq. (16)) and $\hat{H}_{\text{hop}}$ (Eq. (20)) from Sec. III:

[TABLE]

$\tilde{H}_{\text{global}}$ is a model of XY and XZ–plane toric codes $\tilde{H}_{\text{toric}}$ (Eq. (6)) which are coupled together by $\tilde{H}_{t}$ (Eq. (9)) and $\tilde{H}_{h}$ (Eq. (27)), with an infinitesimal perturbation $E^{\text{YZ}}_{\square}\,[\tilde{\sigma},\tilde{\tau}]^{\square_{\text{YZ}}}_{p}$ (Fig. 9) to break accidental degeneracies (discussed in Sec. III.2, IV.1, and IV.2). $\tilde{H}_{\text{toric}}$ and $\tilde{H}_{t}$ were discussed in Sec. III.1. $\tilde{H}_{h}$ and $E^{\text{YZ}}_{\square}\,[\tilde{\sigma},\tilde{\tau}]^{\square_{\text{YZ}}}_{p}$ were discussed in Sec. III.2 in the X-cube basis (Eq. (14)); here we use the decoupled toric code basis (Eq. (5)). $\tilde{H}_{\text{global}}$ is very similar to $\tilde{H}^{\prime}$ (Eq. (28)); the only difference is that $\tilde{H}_{\text{global}}$ replaces $\tilde{H}_{{h^{\prime}}}$ (Eq. (30)) with the $E^{\text{YZ}}_{\square}$ term. (Recall that the choice of $\tilde{H}_{{h^{\prime}}}$ in $\tilde{H}^{\prime}$ was not important.) $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$, $[\tilde{\sigma}]^{+_{\text{XY}}}_{\iota}$ (Fig. 4), and $\tilde{\sigma}^{x}_{\ell}$ are defined on the XY plane toric codes, while $\tilde{\tau}$ operators are defined on XZ plane toric codes. The second line in Eq. (43) is just a short-hand representation of the included terms, and is only shown for convenience. A phase diagram is given in Fig. 10. We will work in a restricted Hilbert space where

[TABLE]

The dual theory looks like a decoupled combination of the two dual theories discussed in Sec. III. Thus, the dual theory consists of a stack of decoupled YZ-plane square lattices with Pauli operators $\eta^{\alpha}_{\hat{\iota}}$ at vertices $\hat{\iota}$ centered on x-axis links $\ell$ of the original cubic lattice, and a grid of weakly coupled Ising chains along the x-axis with Pauli operators $\mu^{\alpha}_{\iota}$ on the same lattice sites $\iota$ as the original cubic lattice. $\eta^{\alpha}_{\hat{\iota}}$ and $\mu^{\alpha}_{\iota}$ will be completely decoupled. The duality is given by

[TABLE]

The mapping is essentially a combination of the duality mappings in Sec. III: Eq. (18) and (23). Both the original theory (modulo the Hilbert space restriction (Eq. (44))) and the dual theory have two qubits per site. The dual Hamiltonian is simply the decoupled sum of the two dual Hamiltonians in Sec. III (Eq. (19) and (24)):

[TABLE]

In Sec. III we described two phase transitions out of the X-cube order phase. $\tilde{H}_{\text{global}}^{\text{dual}}$ describes these transitions in the same way. We now explain how this duality also describes two transitions out of 3+1D $Z_{2}$ topological order.

Let us start in the XY & XZ plane $Z_{2}$ topological order phase with small $t$ and $h$, and then increase $h$ to transition into the 3+1D $Z_{2}$ topological order phase (Fig. 10). In the original model $\tilde{H}_{\text{global}}$ (Eq. (43)), the $h$ coupling hops a composite excitation of XY and XZ–plane toric code charge excitations ($[\tilde{\sigma}]^{+_{\text{XY}}}_{\iota}$ and $[\tilde{\tau}]^{+_{\text{XZ}}}_{\iota}$, see Fig. 1 (c)). In Sec. IV.1 we used degenerate perturbation theory to show that for large hopping strength $h$, this composite charge excitation condenses, and 3+1D $Z_{2}$ topological order results. We now give a physical explanation. After the composite charge condenses, the XY and XZ–plane 2+1D toric code charge excitations are no longer bound to the XY and XZ planes (respectively) since an XY plane charge can split into a XZ plane charge and a composite charge excitation, which means that the XY and XZ–plane charges can freely transform into the other. A XY or XZ–plane charge excitation thus becomes a 3+1D toric code charge excitation after the composite charge condenses. Once the charges are no longer bound to planes, the XY and XZ–plane 2+1D toric code flux excitations are confined since they no longer have well defined braiding statistics with the charge excitations. That is, a pair of XY or XZ–plane flux excitations will incur an energy cost proportional the their distance due to two strings of YZ-plane toric code flux terms $E^{\text{YZ}}_{\square}\,[\tilde{\sigma},\tilde{\tau}]^{\square_{\text{YZ}}}_{p}$ which will also be excited. Thus, the XY and XZ–plane 2+1D toric code flux excitations become 3+1D toric code flux loop excitations when the composite charge condenses.

The dual theory is described by $\tilde{H}_{\text{global}}^{\text{dual}}$ (Eq. (46)). This transition is driven by increasing $h$ (see Fig. 10) while keeping $t$ constant. Since $\tilde{H}^{\text{dual}}$ does not depend on $h$, and is completely decoupled from $\hat{H}_{\text{hop}}^{\text{dual}}$, the phase transition is just described by $\hat{H}_{\text{hop}}^{\text{dual}}$: i.e. decoupled 1+1D Ising chains. The composite charge excitation is dual to an excitation of the 1+1D Ising paramagnet (Fig. 1 (c)). Just like in Sec. III.2, where a similar duality is discussed, $E^{\text{YZ}}_{\square}\,[\tilde{\sigma},\tilde{\tau}]^{\square_{\text{YZ}}}_{p}$ appears to be marginally relevant under RG (see Appendix C) at the phase transition, and the phase transition is likely first order. The extended honeycomb-based model (Eq. (25)) also exhibits this phase transition since $\tilde{H}_{\text{global}}$ (Eq. (43)) is the same as the effective Hamiltonian $\tilde{H}^{\prime}$ (Eq. (28)) of the honeycomb-based model in the $\max(K_{x},K_{y},t,h)\ll K_{z}$ limit (ignoring subleading terms), with the unimportant difference that $\tilde{H}_{\text{global}}$ includes $E^{\text{YZ}}_{\square}\,[\tilde{\sigma},\tilde{\tau}]^{\square_{\text{YZ}}}_{p}$ as a perturbation, while $\tilde{H}^{\prime}$ instead includes $\tilde{H}_{{h^{\prime}}}$ (Eq. (30)).

Let us now consider the 3+1D $Z_{2}$ topological order phase with small $t$ and large $h$, and then increase $t$ to transition into the YZ-plane $Z_{2}$ topological order phase (Fig. 10). In the original model $\tilde{H}_{\text{global}}$ (Eq. (43)), the $t$ coupling creates a YZ-plane loop of four 3+1D toric code fluxes (Fig. 1 (b)). In Sec. IV.2 we used degenerate perturbation theory to show that for strong coupling $t$, the composite flux loops condense, and YZ-plane topological order results. This occurs because when the YZ-plane flux loops condense, the 3+1D toric code charge excitations can no longer move in the x direction, since this would incur a nontrivial braiding statistic with the condensed YZ-plane flux loops. Now consider an open flux loop excitation of 3+1D toric code that is long in the y and z directions, but very short in the x direction. Before the YZ-plane flux loops condense, this excitation has a large energy proportional to its length. However, after the YZ-plane flux loops condense, its bulk is indistinguishable to part of a YZ-plane flux loop, and therefore only a finite energy cost is incurred. Thus, the YZ-plane flux loop condensation allows 2+1D toric code flux excitations (at the end of the open flux loop) to become deconfined along the YZ plane.

The dual theory is described by decoupled YZ-plane 2+1D Ising models $\tilde{H}^{\text{dual}}$ (Eq. (46)) since $\hat{H}_{\text{hop}}^{\text{dual}}$ does not depend on $t$. The YZ-plane flux loop excitation of $\tilde{H}_{\text{global}}$ is dual to a domain wall in a YZ-plane 2+1D Ising model (Fig. 1 (b)). Just like in $Sec.\,\ref{sec:duality}$, where a similar duality is discussed, this phase transition is also likely first order. The extended honeycomb-based model (Eq. (25)) also exhibits this phase transition since $\tilde{H}_{\text{global}}$ (Eq. (43)) is the same as the effective Hamiltonian $\tilde{H}^{\prime}$ (Eq. (28)) of the honeycomb-based model in the $\max(K_{x},K_{y},t,h)\ll K_{z}$ limit (ignoring subleading terms). This $\max(K_{x},K_{y},t,h)\ll K_{z}$ limit describes the phase transition out of 3+1D topological order across the vertical red line in 11. However, this $\max(K_{x},K_{y},t,h)\ll K_{z}$ limit does not correspond to the top-left corner of Fig. 8 for which $\max(K_{x},K_{y},t)\ll\min(K_{z},h)$, which instead corresponds to the diagonal red on the top-left of 11. To describe the $\max(K_{x},K_{y},t)\ll\min(K_{z},h)$ limit, we could apply a similar duality to the following Hamiltonian, which describes the extended honeycomb-based model (Eq. (25)) when $\max(K_{x},K_{y},t)\ll\min(K_{z},h)$:

[TABLE]

where $\hat{H}_{\text{3D toric}}$ is given in Eq. (36).

IV.4 Connections to Previous Coupled Layer Constructions

$\tilde{H}_{\text{toric}}+\tilde{H}_{t}+\tilde{H}_{h}$ (Eq. (6), (9), and (30)) is very similar to the coupled layer constructions considered in Ref. Vijay (2017) and Ma et al. (2017). The primary difference is that we have only XY and XZ–plane toric codes, while they also had YZ planes. This difference was motivated in our work by the desire to have the simplest possible lattice model. Although the lack of YZ planes in $\tilde{H}_{\text{toric}}$ adds some anisotropy to the resulting X-cube Hamiltonian (Eq. (15)) and 3+1D toric code Hamiltonian (Eq. (36)), these details do not affect the universality class of these phases. However, when both $t$ and $h$ are large, the lack of YZ-plane toric codes in our model results in radically different physics. Specifically, in Sec. IV.2 we showed that XY and XZ–plane toric codes with large $t$ and $h$ coupling leads to a phase described by decoupled YZ-plane toric codes. However, with XY, XZ, and YZ planes, Ref. Vijay (2017) showed that large $t$ and $h$ coupling drives the system into a trivial phase with no topological order. It is interesting that large $t$ and $h$ coupling seems to have a $Z_{2}$ action on the presence of topological order on YZ planes.

V Conclusion

In this work we strive to bring the recently discovered fracton topological order closer to an experimental realization. We do this by presenting a simple model with X-cube fracton topological order which only requires nearest-neighbor 2-spin interactions. This is a significant improvement over the original X-cube model with 12-spin interactions Vijay et al. (2016), or previous coupled layer constructions with 4-spin interactions. Unfortunately however, it is likely that much future work will be necessary to realize fracton order in an experiment.

Perhaps there exists a material for which our model could be a useful description. If this is not the case, we hope that this work could be useful inspiration for deriving more realistic models of fracton order with possible experimental realizations. Additionally, synthetic quantum matter, such as cold atom experiments, provide another possible route to experimental realization. Perhaps previous theoretical work Tanamoto et al. (2015); Liu and Xue (2013) on synthetic quantum matter realizations of Kitaev’s honeycomb model Kitaev (2006) could be generalized to apply to our model or to similar models. It would also be interesting to discover realistic models of other kinds of fracton order. In particular, it would be very interesting to realize fractal-like (type II in Ref. Vijay et al. (2016)) fracton topological order (e.g. Haah’s code Haah (2011)) via a coupled layer construction or via a two-spin Hamiltonian.

In this work, we found two first order phase transitions out of the X-cube phase. Ref. Vijay (2017) also showed that their isotropic coupled layer construction admits dualities describing two first order phase transition out of X-cube fracton topological order. Additionally, Ref. Vijay et al. (2016) gave dualities describing phase transitions out of an entire class of fracton phases; however, the nature of these phase transitions is currently unknown. It remains an open and interesting question if it is possible to have a continuous phase transition out of fracton topological order.

In Fig. 8 and especially in Fig. 3, there are regions of the phase diagrams which could not be described by degenerate perturbation theory or the dualities that we found. In particular, it is not clear what occurs when $K_{x}$ and $K_{y}$ are increased above the $\max(K_{x},K_{y})\ll K_{z}$ limit in Fig. 3. A gapless version of X-cube fracton order would be an interesting possibility; however, such a phase is likely to be unstable to a gapped phase Vijay (2017); Xu and Wu (2008). It is also natural to ask what occurs near the center of Fig. 8.

It is also interesting to note that the X-cube operators that we generated using degenerate perturbation theory ($[\sigma,\tau]^{\text{\mancube}}_{p}$ (Fig. 6), and $[\sigma]^{+_{\text{XY}}}_{\iota}$ and $[\tau]^{+_{\text{XZ}}}_{\iota}$ (Fig. 4)) all commute with the original Hamiltonian (Eq. (1)). The original Hamiltonian has 8 qubit degrees of freedom for each unit cell; and we have found 3 qubit integrals of motion. In Kitaev’s exactly solvable honeycomb model Kitaev (2006), every unit cell is comprised of 2 qubit degrees of freedom which are split into a single integral of motion given by a hexagon/vison operator (Fig. 4), and two Majorana fermions with only quadratic interactions for a given static configuration of hexagon/vison excitations. Perhaps our model, or a similar model, admits a similar exact solution, possibly with connections to the generalized parton constructions in Ref. Hsieh and Halász (2017).

Acknowledgements.

Acknowledgements — We thank Leon Balents for pointing out the relevance of the layer couplings discussed in Appendix C. This work was supported by the NSERC of Canada and the Center for Quantum Materials at the University of Toronto. Computations were performed on the gpc supercomputer at the SciNet HPC Consortium. SciNet is funded by: the Canada Foundation for Innovation under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund - Research Excellence; and the University of Toronto. Loken et al. (2010)

Appendix A Large $K_{z}$ Degenerate Perturbation Theory

In this section, we will apply degenerate perturbation theory (DPT) (Appendix B) to the extended model $H^{\prime}$ (Eq. (25)) and produce a phase diagram (Fig. 11) in the limit when $K_{z}\gg\max(K_{x},K_{y})$. To reduce verbosity, throughout this section we will often refer to operators by their coefficient. For example, when we say that $E_{\square}$ and $t$ anticommute, we mean that an operator with coefficient $E_{\square}$ anticommutes with an operator with coefficient $t$. We will perform DPT by splitting the Hamiltonian in Eq. (25) into $H^{\prime}=H_{0}+H_{1}$ where $|H_{0}|\gg|H_{1}|$ (where any operator norm can be chosen). $K_{z}$ will always be included in $H_{0}$, and $K_{x}$ and $K_{y}$ will always be in $H_{1}$. Therefore, $t$ can only be in $H_{0}$ if $t\gg\max(K_{x},K_{y})$; and $t$ can only be in $H_{1}$ if $t\ll K_{z}$. Similarly, $h$ can only be in $H_{0}$ if $h\gg\max(K_{x},K_{y})$; and $h$ can only be in $H_{1}$ if $h\ll K_{z}$. These (overlapping) conditions partition (via the dotted gray lines) the phase diagram (Fig. 11) into four overlapping quadrants.

DPT will generate many different terms. But how do we determine the resulting phase? Instead of trying to answer this question direction, we will instead consider one target phase at a time (e.g. X-cube fracton order), and then attempt to calculate a sufficient condition for that phase (e.g. $\max(K_{x},K_{y})\ll t\ll K_{z}$). In our model, doing this is simplified because all terms in the Hamiltonian that we consider are products of Pauli operators, which we will refer to as Pauli strings. Furthermore, all of the phases that we’re interested in are gapped and can be described by a stabilizer Hamiltonian. By stabilizer Hamiltonian, we will mean a Hamiltonian that is a sum of Pauli strings which all commute with each other. The ground state of the stabilizer Hamiltonian will also be the ground state of an extensive number of Pauli string terms in the stabilizer Hamiltonian; we will refer to these Pauli strings as stabilizers. Examples of stabilizer Hamiltonians are: $H=\sum_{i}X_{i}$ which stabilizes a paramagnet, $H=\sum_{\langle ij\rangle}Z_{i}Z_{j}$ which stabilizes a ferromagnet, or the toric code $H=\sum_{\square}\prod_{i\in\square}Z_{i}+\sum_{+}\prod_{i\in+}X_{i}$ which stabilizes $Z_{2}$ topological order. In certain limits, such as $\max(K_{x},K_{y})\ll t\ll K_{z}$ in Eq. (7), we will find that (multiple applications of) DPT produces a Hamiltonian that is a stabilizer Hamiltonian. If we obtain the right stabilizers, then we know that the limits considered are a sufficient condition for the target phase.

As a generic example, suppose we want to find sufficient conditions for a target phase with stabilizers: $A_{i}$, $B_{i}$, $C_{i}$, and $D_{i}$. If at least one of these operators appears in $H$, then we can apply DPT with at least one of them appearing in $H_{0}$. Suppose only $A_{i}$ appears in $H$ with energy $E_{A}$. Then we must choose $H_{0}=-\sum_{i}E_{A}A_{i}$. In order for DPT to converge, we will have to assume that $E_{A}$ is much larger than all energies in $H_{1}$. This will thus be part of the sufficient condition for the target phase. If we assume this, and generate an effective Hamiltonian $H_{\text{eff}}$, then $A_{i}$ will be a stabilizer (after more applications of DPT) of the ground state since $H_{\text{eff}}$ resides in the degenerate ground state manifold of $H_{0}$. Now suppose that $\sum_{i}E_{B}B_{i}$ is a term in $H_{\text{eff}}$, but neither $C_{i}$ nor $D_{i}$ appear. Then we will want to apply DPT to $H_{\text{eff}}$ with $H_{0}^{\prime}=-\sum_{i}E_{B}B_{i}$ and $H_{1}^{\prime}=H_{\text{eff}}-H_{0}^{\prime}$ to generate a new effective Hamiltonian $H_{\text{eff}}$ in hope of eventually generating $C_{i}$ and $D_{i}$. In order for DPT to converge, we will have to assume that $E_{B}$ is larger than all other terms in $H_{\text{eff}}$. In particular, we will be most worried about terms which anticommute with $B_{i}$, since these terms will create excitations of $B_{i}$ and will therefore try to prevent $B_{i}$ from being a stabilizer. Suppose the new effective Hamiltonian $H_{\text{eff}}^{\prime}$ does include $C_{i}$ and $D_{i}$ with energies $E_{C}$ and $E_{D}$, respectively. Then if $E_{C}$ and $E_{D}$ are sufficiently large, we can apply DPT once more with $H_{0}^{\prime\prime}=-\sum_{i}(E_{C}C_{i}+E_{D}D_{i})$, and we will be left with a stabilizer Hamiltonian of the desired phase. Note that terms like $\sum_{\langle ij\rangle}E_{CC}C_{i}C_{j}$ may also be generated. However, as long as $E_{CC}\ll E_{C}$, these terms don’t affect the ground state, and so they will be ignored since $C_{i}C_{j}$ is not an independent stabilizer from $C_{i}$ (since it’s the product of $C_{i}$ and $C_{j}$). Thus, a sufficient condition for the target phase would require that $E_{A}$, $E_{B}$, $E_{C}$, and $E_{D}$ are sufficiently large.

In Fig. 11 we show the results from this approach. The sufficient conditions that we generate for each phase are not necessary conditions. Thus, we can not fill in the entire phase diagram (as seen in Fig. 8). However, Fig. 11 is drawn in a subspace of couplings ($K_{\mu},t,h$) where we can identify the phase at almost every point.

Below we list the energies (calculated in Appendix A.5) for the important terms generated by degenerate perturbation theory (DPT) that are required for drawing the phase diagram (Fig. 11).

[TABLE]

$E_{\square}$ (introduced in Eq. (4)) corresponds to a toric code flux operator (Fig. 4) on an XY or XZ plane. $E_{+}$ (Eq. (4)) corresponds to a 2+1D toric code charge operator (Fig. 4) on an XY or XZ plane. The “$\sim$” in the expressions means that these expressions are only asymptotically correct; e.g. factors of two are ignored which allows simplifications such as $K_{z}+t\sim\max(K_{z},t)$.

Energies for other terms generated by DPT will also be mentioned. These energies can be difficult to calculate in some cases due to the high order of DPT required to calculate them; but fortunately their values will not be important in this work. These energies are:

•

$E^{\text{YZ}}_{\square}$ (in Eq. (33) and (39)), corresponding to a toric code flux operator (Fig. 9) on a YZ plane

•

$E_{\text{\mancube}}$ (in Eq. (12)), corresponding to an X-cube fracton number operator (Fig. 6)

•

$E_{*}$ (in Eq. (33)), corresponding to a 3+1D toric code charge operator (Fig. 9)

•

$E^{\text{YZ}}_{+}$ (in Eq. (39)), corresponding to the same operator as $E_{*}$, but we instead view it as a YZ-plane toric code charge operator

We now apply DPT to $H^{\prime}$ (Eq. (25)) in various limits. $H_{{h^{\prime}}}$ (Eq. (27)) will be considered as an example of a generic perturbation that can be added to break unphysical degeneracy which results in some cases; it’s coefficient ${h^{\prime}}$ will always be assumed to be infinitesimally small. Then we will calculate $E_{\square}$ and $E_{+}$ (Eq. (48)) in Appendix A.5.

A.1 XY&XZ plane 2+1D $Z_{2}$ Topological Order

This phase occurs when $K_{z}$, $E_{\square}$, and $E_{+}$ are the stabilizers. Let us determine sufficient conditions for this. We’ll need to do DPT in order to produce $E_{\square}$ and $E_{+}$. We’ll do DPT with $H_{0}=-K_{z}\,(\cdots)$. To ensure that DPT converges, we’ll need to assume $K_{z}\gg\max(K_{x},K_{y},t,h)$. $E_{\square}$ and $E_{+}$ are the highest energy operators generated by DPT. But $E_{\square}$ anticommutes with $t$, and can therefore only be a stabilizer if $E_{\square}\gg t$. And $E_{+}$ anticommutes with $h$, and can therefore only be a stabilizer if $E_{+}\gg h$. Therefore, if

[TABLE]

then another application of DPT with $H_{0}^{\prime}=-K_{z}\,(\cdots)-E_{\square}\,(\cdots)-E_{+}\,(\cdots)$ will produce a stabilizer Hamiltonian with 2+1D $Z_{2}$ topological order along the original XY and XZ planes. The above conditions are already sufficient to ensure that DPT converges. This is the decoupled phase.

A.2 X-cube Fracton Topological Order

This phase occurs when $K_{z}$, $t$, $E_{\text{\mancube}}$, and $E_{+}$ are stabilizers. Let us determine sufficient conditions for this. We’ll need to do DPT in order to produce $E_{\text{\mancube}}$ and $E_{+}$.

Let us first try DPT with $H_{0}=-K_{z}\,(\cdots)-t\,(\cdots)$. To ensure that DPT converges, we’ll need to assume $\min(K_{z},t)\gg\max(K_{x},K_{y},h)$. $E_{+}$ is the highest energy operator generated by DPT. But $E_{+}$ anticommutes with $h$, and can therefore only be a stabilizer if $E_{+}\gg h$. $E_{\text{\mancube}}$ is the next highest energy operator generated by DPT that’s independent from and commutes with the previous (candidate) stabilizers: $K_{z}$, $t$, and $E_{+}$. (The product of two $E_{+}$ terms may be generated with a higher energy, but such terms will not affect the ground state since their energy is much smaller than $E_{+}$.) $E_{\text{\mancube}}$ can therefore be a stabilizer. In fact, if $K_{z}$, $t$, and $E_{+}$ are stabilizers, then $E_{\text{\mancube}}$ will have to be a stabilizer if we require the ground state to be the ground state of a stabilizer Hamiltonian which is stable to arbitrary perturbations. There simply is no other Pauli string that commutes with $K_{z}$, $t$, and $E_{+}$, and isn’t a product of $K_{z}$, $t$, $E_{+}$, and $E_{\text{\mancube}}$. Therefore, if $\min(K_{z},t)\gg\max(K_{x},K_{y},h)$ and $E_{+}\gg h$, then X-cube fracton topological order will be stabilized. That is, applying DPT again with $H_{0}^{\prime}=-K_{z}\,(\cdots)-t\,(\cdots)-E_{\text{\mancube}}\,(\cdots)-E_{+}\,(\cdots)$ will result in a stabilizer Hamiltonian with X-cube fracton order; the above assumptions are already sufficient to guarantee that DPT converges.

But what if it isn’t the case that $t\gg\max(K_{x},K_{y},h)$? To cover this case, we’ll also try DPT with $H_{0}=-K_{z}\,(\cdots)$. To ensure that DPT converges, we’ll need to assume $K_{z}\gg\max(K_{x},K_{y},t,h)$. (This assumption will overlap with the $\min(K_{z},t)\gg\max(K_{x},K_{y},h)$ assumption in the previous paragraph, but the results in this overlapping region are the same.) $E_{\square}$ and $E_{+}$ are the highest energy operators generated by DPT. $t$ anticommutes with $E_{\square}$, and can therefore only be a stabilizer if $t\gg E_{\square}$. We can do DPT again with $H_{0}^{\prime}=-t\,(\cdots)$ and $H_{1}^{\prime}=-E_{\square}\,(\cdots)-E_{+}\,(\cdots)-h\,(\cdots)+\cdots$ where the last “$\cdots$” denotes unimportant subleading terms. $t\gg E_{\square}$ (note $E_{\square}\sim E_{+}$) is already sufficient to ensure that DPT converges. $E_{+}$ anticommutes with $h$, and can therefore only be a stabilizer if $E_{+}\gg h$. $E_{\text{\mancube}}$ is the highest energy operator generated by DPT that’s independent from and commutes with the previous (candidate) stabilizers: $K_{z}$, $t$ and $E_{+}$. $E_{\text{\mancube}}$ can therefore be a stabilizer. Therefore, if $K_{z}\gg\max(K_{x},K_{y},t,h)$ and $t\gg E_{\square}$ and $E_{+}\gg h$, then X-cube fracton topological order will be stabilized.

In summary, a sufficient condition for X-cube fracton topological order is

[TABLE]

A.3 3+1D $Z_{2}$ Topological Order

This phase occurs when $K_{z}$, $h$, $E_{\square}$, $E^{\text{YZ}}_{\square}$, and $E_{*}$ are the stabilizers. Let us determine sufficient conditions for this. We’ll need to do DPT in order to produce $E_{\square}$, $E^{\text{YZ}}_{\square}$, and $E_{*}$.

Let us first try DPT with $H_{0}=-K_{z}\,(\cdots)-h\,(\cdots)$. To ensure that DPT converges, we’ll need to assume $\min(K_{z},h)\gg\max(K_{x},K_{y},t)$. $E_{\square}$ is the highest energy operator generated by DPT. But $E_{\square}$ anticommutes with $t$, and can therefore only be a stabilizer if $E_{\square}\gg t$. $E_{\text{\mancube}}$ and $E_{*}$ are the next highest energy operators generated by DPT that are independent from and commute with the previous (candidate) stabilizers: $K_{z}$, $h$, and $E_{\square}$. $E_{\text{\mancube}}$ and $E_{*}$ can therefore be stabilizers since they commute. (If they anticommuted, then they couldn’t both be stabilizers.) The current stabilizers ($K_{z}$, $h$, $E_{\square}$, $E_{*}$, and $E_{\text{\mancube}}$) describe a phase corresponding to 3+1D toric code, but without YZ-plane plaquette/flux operators. This phase is unstable to perturbations and has unphysical degeneracy $\sim 2^{L^{2}}$ on a torus due to the $L^{2}$ missing $E^{\text{YZ}}_{\square}$ operators on any single YZ plane. However, with generic perturbations (for which we choose $H_{{h^{\prime}}}$ (Eq. (27)) as an example), we can do DPT again with $H_{0}^{\prime}=-E_{\square}\,(\cdots)-E_{*}\,(\cdots)-E_{\text{\mancube}}\,(\cdots)$ and $H_{1}^{\prime}=H_{{h^{\prime}}}+\cdots$ where the last “$\cdots$” denotes unimportant subleading terms. $E^{\text{YZ}}_{\square}$ is the highest energy operator generated by DPT. $E^{\text{YZ}}_{\square}$ can therefore be a stabilizer (which makes $E_{\text{\mancube}}$ a redundant stabilizer since it’s the product of two $E^{\text{YZ}}_{\square}$ operators, modulo $K_{z}$ and $h$ operators). Therefore, if $\min(K_{z},h)\gg\max(K_{x},K_{y},t)$ and $E_{\square}\gg t$, then 3+1D $Z_{2}$ topological order will be stabilized.

But what if it isn’t the case that $h\gg\max(K_{x},K_{y},t)$? To cover this case, we’ll also try DPT with $H_{0}=-K_{z}\,(\cdots)$. To ensure that DPT converges, we’ll need to assume $K_{z}\gg\max(K_{x},K_{y},t,h)$. (This assumption will overlap with the $\min(K_{z},h)\gg\max(K_{x},K_{y},t)$ assumption in the previous paragraph, but the results in this overlapping region are the same.) $E_{\square}$ and $E_{+}$ are the highest energy operators generated by DPT. But $h$ anticommutes with $E_{+}$, and can therefore only be a stabilizer if $h\gg E_{+}$. We can do DPT again with $H_{0}^{\prime}=-h\,(\cdots)$ and $H_{1}^{\prime}=-E_{\square}\,(\cdots)-E_{+}\,(\cdots)-t\,(\cdots)+\cdots$. $h\gg E_{+}$ is already sufficient to ensure that DPT converges. $E_{\square}$ anticommutes with $t$, and can therefore only be a stabilizer if $E_{\square}\gg t$. $E_{\text{\mancube}}$ and $E_{*}$ are the highest energy operators generated by DPT that are independent from and commute with the previous (candidate) stabilizers: $K_{z}$, $h$ and $E_{\square}$. $E_{\text{\mancube}}$ and $E_{*}$ can therefore be stabilizers since they commute. We can do DPT again with $H_{0}^{\prime\prime}=-h\,(\cdots)-E_{\square}\,(\cdots)-E_{*}\,(\cdots)-E_{\text{\mancube}}\,(\cdots)$ and $H_{1}^{\prime\prime}=H_{{h^{\prime}}}+\cdots$. $E^{\text{YZ}}_{\square}$ is the highest energy operator generated by DPT and can therefore be a stabilizer. Therefore, if $K_{z}\gg\max(K_{x},K_{y},t,h)$ and $h\gg E_{+}$ and $E_{\square}\gg t$, then 3+1D topological order will be stabilized.

In summary, a sufficient condition for 3+1D $Z_{2}$ topological order is

[TABLE]

when generic (but arbitrarily small) perturbations (such as $H^{\prime}$ (Eq. (25))) are included in the Hamiltonian $H^{\prime}$ (Eq. (25)).

A.4 YZ-plane 2+1D $Z_{2}$ topological order

This phase occurs when $K_{z}$, $h$, $t$, $E^{\text{YZ}}_{\square}$, and $E^{\text{YZ}}_{+}$ are the stabilizers. Let us determine sufficient conditions for this. We’ll need to do DPT in order to produce $E^{\text{YZ}}_{\square}$ and $E^{\text{YZ}}_{+}$.

Let us first try DPT with $H_{0}=-K_{z}\,(\cdots)-t\,(\cdots)+h\,(\cdots)$. To ensure that DPT converges, we’ll need to assume $\min(K_{z},t,h)\gg\max(K_{x},K_{y})$. $E^{\text{YZ}}_{+}$ and $E_{\text{\mancube}}$ are the highest energy operators generated by DPT. $E^{\text{YZ}}_{+}$ and $E_{\text{\mancube}}$ can therefore be stabilizers since they commute. The current stabilizers ($K_{z}$, $h$, $E^{\text{YZ}}_{+}$, and $E_{\text{\mancube}}$) describe a phase that is unstable to perturbations and has unphysical degeneracy $\sim 2^{L^{2}}$ on a torus due to the $L^{2}$ missing $E^{\text{YZ}}_{\square}$ operators on any single YZ plane. However, with generic perturbations (for which we choose $H_{{h^{\prime}}}$ (Eq. (27)) as an example), we can do DPT again with $H_{0}^{\prime}=-E^{\text{YZ}}_{+}\,(\cdots)-E_{\text{\mancube}}\,(\cdots)$ and $H_{1}^{\prime}=H_{{h^{\prime}}}+\cdots$ where the last “$\cdots$” denotes unimportant subleading terms. $E^{\text{YZ}}_{\square}$ is the highest energy operator generated by DPT and can therefore be a stabilizer. Therefore, if $\min(K_{z},t,h)\gg\max(K_{x},K_{y})$, then the desired phase will be stabilized.

Let us now try DPT with $H_{0}=-K_{z}\,(\cdots)-t\,(\cdots)$. To ensure that DPT converges, we’ll need to assume $\min(K_{z},t)\gg\max(K_{x},K_{y},h)$. $E_{+}$ is the highest energy operator generated by DPT. But $h$ anticommutes with $E_{+}$, and can therefore only be a stabilizer if $h\gg E_{+}$. $E_{+}$ and $E_{\text{\mancube}}$ are the next highest energy operators generated by DPT that are independent from and commute with the previous (candidate) stabilizers: $K_{z}$, $t$, and $h$. $E_{+}$ and $E_{\text{\mancube}}$ can therefore be stabilizers since they commute. We can do DPT again with $H_{0}^{\prime}=-h\,(\cdots)-E_{+}\,(\cdots)-E_{\text{\mancube}}\,(\cdots)$ and $H_{1}^{\prime}=H_{{h^{\prime}}}+\cdots$. $E^{\text{YZ}}_{\square}$ is the highest energy operator generated by DPT and can therefore be a stabilizer. Therefore, if $\min(K_{z},t)\gg\max(K_{x},K_{y},h)$ and $h\gg E_{+}$, then the desired phase will be stabilized.

Let us now try DPT with $H_{0}=-K_{z}\,(\cdots)-h\,(\cdots)$. To ensure that DPT converges, we’ll need to assume $\min(K_{z},h)\gg\max(K_{x},K_{y},t)$. $E_{\square}$ is the highest energy operator generated by DPT. But $t$ anticommutes with $E_{\square}$, and can therefore only be a stabilizer if $t\gg E_{\square}$. $E^{\text{YZ}}_{+}$ and $E_{\text{\mancube}}$ are the next highest energy operators generated by DPT that are independent from and commute with the previous (candidate) stabilizers: $K_{z}$, $t$, and $h$. $E^{\text{YZ}}_{+}$ and $E_{\text{\mancube}}$ can therefore be stabilizers since they commute. We can do DPT again with $H_{0}^{\prime}=-t\,(\cdots)-E^{\text{YZ}}_{+}\,(\cdots)-E_{\text{\mancube}}\,(\cdots)$ and $H_{1}^{\prime}=H_{{h^{\prime}}}+\cdots$. $E^{\text{YZ}}_{\square}$ is the highest energy operator generated by DPT and can therefore be a stabilizer. Therefore, if $\min(K_{z},h)\gg\max(K_{x},K_{y},t)$ and $t\gg E_{\square}$, then the desired phase will be stabilized.

Let us now try DPT with $H_{0}=-K_{z}\,(\cdots)$. To ensure that DPT converges, we’ll need to assume $K_{z}\gg\max(K_{x},K_{y},t,h)$. $E_{\square}$ and $E_{+}$ are the highest energy operators generated by DPT. But $t$ anticommutes with $E_{\square}$, and can therefore only be a stabilizer if $t\gg E_{\square}$. And $h$ anticommutes with $E_{+}$ , and can therefore only be a stabilizer if $h\gg E_{+}$. We can do DPT again with $H_{0}^{\prime}=-t\,(\cdots)-h\,(\cdots)$ and $H_{1}^{\prime}=-E_{\square}\,(\cdots)-E_{+}\,(\cdots)+\cdots$ where the last “$\cdots$” denotes unimportant subleading terms. The previous assumptions are sufficient to ensure that DPT converges. $E_{\text{\mancube}}$ and $E^{\text{YZ}}_{+}$ are the highest energy operators generated by DPT that are independent from and commute with the previous (candidate) stabilizers: $K_{z}$, $t$, and $h$. $E_{\text{\mancube}}$ and $E^{\text{YZ}}_{+}$ can therefore be stabilizers since they commute. We can do DPT again with $H_{0}^{\prime\prime}=-t\,(\cdots)-h\,(\cdots)-E^{\text{YZ}}_{+}\,(\cdots)-E_{\text{\mancube}}\,(\cdots)$ and $H_{1}^{\prime\prime}=H_{{h^{\prime}}}+\cdots$. $E^{\text{YZ}}_{\square}$ is the highest energy operator generated by DPT and can therefore be a stabilizer. Therefore, if $K_{z}\gg\max(K_{x},K_{y},t,h)$ and $t\gg E_{\square}$ and $h\gg E_{+}$, then the desired phase will be stabilized.

In summary, a sufficient condition for 2+1D $Z_{2}$ topological order on YZ planes is

[TABLE]

when generic (but arbitrarily small) perturbations (such as $H^{\prime}$ (Eq. (25))) are included in the Hamiltonian $H^{\prime}$ (Eq. (25)). Physically, when these perturbations are small, $E_{\text{\mancube}}\gg E^{\text{YZ}}_{\square}$ and this implies that in the excited states, flux excitations $E^{\text{YZ}}_{\square}$ will want to be bound to a partner flux excitation on a neighboring layer. However, this does not change the resulting phase of the ground state.

A.5 $E_{\square}$ and $E_{+}$ Calculation

In this section, we’ll calculate $E_{\square}$ and $E_{+}$ using degenerate perturbation theory (Appendix B). $E_{\square}$ and $E_{+}$ appear in the effective Hamiltonian $H^{\text{eff}}$ at fourth order. and the only term at this order that contributes is $\mathcal{P}(H_{1}\mathcal{D})^{3}H_{1}\mathcal{P}$. We will focus on the $E_{\square}$ or $E_{+}$ term that is generated for a single hexagon with indices labeled as shown below.

Thus, we only need to consider the terms of $H_{1}$ that act only on these six sites. And we can erase all Pauli operators from $H_{0}$ which aren’t on one of these six sites. $E_{\square}$ is centered on plaquettes (Fig. 4) and is calculated assuming $t\ll K_{z}$. Thus, we can take $t=0$ when calculating $E_{\square}$, but we must consider the $h$ term acting on sites 1 and 4, and 3 and 6. $E_{+}$ is centered on vertices (Fig. 4) and is calculated assuming $h\ll K_{z}$. Thus, we can take $h=0$ when calculating $E_{+}$, but we must consider the $t$ term acting on sites 2 and 5. Therefore, we can calculate $E_{\square}$ and $E_{+}$ at the same time by using the $H_{0}$ and $H_{1}$ below and taking either $t$ or $h$ equal to zero.

[TABLE]

The effective Hamiltonian is:

[TABLE]

$E_{\square}$ is given by $E$ in Eq. (54) with $t=0$, which matches Eq. (48) asymptotically. $E_{+}$ is given by $E$ in Eq. (54) with $h=0$, which also matches Eq. (48) asymptotically. Note that we used the fact that $\sigma^{z}_{1}\sigma^{z}_{4}=\sigma^{z}_{3}\sigma^{z}_{6}=1$ in the low energy Hilbert space to get Eq. (53) from the previous line, which matches the form of $[\sigma]^{\square_{\text{XY}}}_{p}$ and $[\sigma]^{+_{\text{XY}}}_{\iota}$ in Fig. 4.

Appendix B Degenerate Perturbation Theory

In this section, we review the Schrieffer-Wolff formulation of degenerate perturbation theory Schrieffer and Wolff (1966); Bravyi et al. (2011b), which we use in the rest of the paper. We begin by splitting a Hamiltonian $H$ into an unperturbed Hamiltonian $H_{0}$ and a perturbation $\lambda H_{1}$ where $\lambda\ll 1$.

[TABLE]

We will focus on the case where $H_{0}$ has a degenerate ground state manifold with energy $E_{0}$. The goal is to derive an effective Hamiltonian $H^{\text{eff}}$ which describes the physics of the ground state manifold in the presence of the perturbation $\lambda H_{1}$.

We define an operator $\mathcal{P}$ which projects onto the ground state manifold:

[TABLE]

and define the effective Hamiltonian $H^{\text{eff}}$ using a canonical transformation on $H$

[TABLE]

The canonical transformation is generated by an antihermitian operator $S$. We also define $X_{n}$ as shown above. With $S$ defined as in Eq. (58), in order for $H^{\text{eff}}$ to capture the desired physics, it is sufficient to require that the effective Hamiltonian commutes with the ground state projection operator $\mathcal{P}$:

[TABLE]

This condition implies that the unitary rotation $e^{S}$ has rotated $H$ into $H^{\text{eff}}$ with a block diagonal form such that the original ground state manifold is completely decoupled from the high energy Hilbert space. A solution to Eq. (59) exists in which $S_{n}$ is defined iteratively in terms of $S_{1},S_{2},\cdots,S_{n-1}$:

[TABLE]

The disadvantage of this solution is that $S$ is a nonlocal operator due to the presence of the nonlocal projection operator $\mathcal{P}$. Thus, when we calculate a series expansion for $H^{\text{eff}}$ (Eq. (63)), it will not be explicitly local. The locality of $H^{\text{eff}}$ will only result from a precise cancellation due to minus signs.

We now define some compact notation:

[TABLE]

Using this notation and iterating Eq. (60) on a computer gives the following expansion for the canonical transformation generator $S$ and the effective Hamiltonian $H^{\text{eff}}$ projected into the low energy Hilbert space:

[TABLE]

[TABLE]

Appendix C Monte Carlo of Dual Ising Theories

In this appendix, we use classical Monte Carlo to provide evidence that the weakly coupled stacks of transverse-field Ising models that describe the dual phase transitions in Fig. 1 are likely to be first order. There are two cases to consider, 1) a weakly coupled 1D stack of 2+1D Ising models (as in Sec. III.1), and 2) a weakly coupled 2D grid of 1+1D Ising models (as in Sec. III.2). A classical Ising model in 3D and 2D will be used to simulate the 2+1D and 1+1D (respectively) zero temperature transverse field Ising model layers.

$I$ and $J$ will be used to denote the 3D or 2D lattice coordinate, while $L$ will denote the layer (of the 1D stack or 2D grid). Thus, we will be interested in the following classical layered Ising-like model (with inverse temperature $\beta=1$):

[TABLE]

$\sum_{L}$ sums over the layers and $\sum_{\langle IJ\rangle}$ sums over nearest-neighbor sites. The $J$ term is a nearest-neighbor ferromagnetic Ising coupling that only couples spins on the same layer, where $J$ is roughly related by $J\sim E_{\square}/t$ in Eq. (19) or $J\sim h/E_{+}$ in Eq. (24). The $K$ term is a 4-spin term which couples together neighbor-neighbor spins or neighboring layers, and roughly corresponds to a product of two $[\tilde{\sigma}]^{\square_{\text{XY}}}_{p}$ or two $[\tilde{\tau}]^{\square_{\text{XZ}}}_{p}$ operators in Sec. III.1, or $J\sim E^{\text{YZ}}_{\square}/E_{+}$ in Eq. (24). It is useful to rewrite $H_{\text{Ising}}$ as in Eq. (65) to reveal that the $K$ term results in a new effective $J$ term ($J_{LIJ}^{\text{eff}}$) which is site and layer dependent. $\sum_{\langle LL^{\prime}\rangle}^{L}$ sums over all layers $L^{\prime}$ which neighbor the layer $L$.

If $K=0$, then the layers are decoupled and every layer independently undergoes an order to disorder transition. If $K$ is small and positive, then we can see from the form of $J_{LIJ}^{\text{eff}}$ that when a layer transitions into its ordered state, the average $J_{LIJ}^{\text{eff}}$ for neighboring layers will increase. Thus, positive $K$ causes the effective $J$ to increase rapidly at the disorder to order transitions when $J$ is increased. It is reasonable that this could lead to a first order transition.

If $K$ is small and negative, then when a layer transitions into its ordered state, the average $J_{LIJ}^{\text{eff}}$ for neighboring layers will decrease. This makes it difficult for the system to order all at once. What happens instead for intermediate $J$ is that only half the layers order such that ordered layers neighbor disordered layers, and visa versa. Thus, $K$ makes the average $J_{LIJ}^{\text{eff}}$ smaller only for the disordered layers. In Sec. C.1, we will briefly discuss what these phases correspond to in the original model.

From the form of Eq. (65) we see that each layer is an Ising model with a site-dependent coupling $J_{LIJ}^{\text{eff}}$. Therefore, we can apply the global Wolff update algorithm Wolff (1989) to this model one layer at a time. Thus, it is easy to generalize Wolff update code to efficiently study this model even when the correlation length along the layers is large. However, since the update is only applied to one layer at a time, the update is effectively local across the layers. Thus, if the correlation length across the layers is large, the simulation will be inefficient. This can occur when the layer coupling $K$ is large; however, since we are interested in weak layer coupling $K$, this isn’t a major issue.

After applying this algorithm Slagle (2017a), we find that all of the transitions appear to be first order. The simulation results are summarized in Fig. 12. The figure shows that at every phase transition, $\xi_{c}/L$ tends to decrease as the system size $L$ increases, where $\xi_{c}$ is the connected correlation length between two spins on the same layer. 111The connected correlation length $\xi_{c}$ was calculated using Eqn. 73 in Sandvik (2011) (At a continuous phase transition, $\xi_{c}/L$ would instead asymptote to a constant.) We also find other signs of discontinuous phase transitions, such as magnetization histograms with two bumps, which indicates the presence of metastable states and hysteresis effects which occur in discontinuous transitions. However, we can not rule out the possibility that smaller $K$ could instead result in a continuous phase transition. Indeed, our smaller $K$ looks second order at the currently simulated lattice sizes; however, we expect that larger lattices would then likely show signs of discontinuous phase transitions.

C.1 Intermediate Phases

From Fig. 12 we see that when the layer coupling $K$ is negative, there is an intermediate phase where the squared layer magnetization $\langle S_{L}^{2}\rangle\sim 1/2$, where $S_{L}\equiv\frac{1}{N_{I}}\sum_{I}\sigma_{LI}$ is the average magnetization on a given layer $L$. As previously mentioned, this phase spontaneously breaks translation symmetry and is characterized by a checkerboard pattern where half of the layers order, while the other half remain disordered. This can be described by the following order parameter

[TABLE]

In Fig. 13 we show that this order parameter $R$ is nonzero only in the intermediate phase.

These intermediate phases of the layered Ising models are dual phases of phases in the original fracton model. We can determine which phase results in the original model by adding a term to the dual Ising model which explicitly breaks the symmetry that is spontaneously broken in the intermediate phase, and then apply the appropriate duality mapping (Eq. (18), 23 or 45) in reverse. To explicitly break the symmetry, we add the following term to $H_{\text{Ising}}$ (Eq. (64)):

[TABLE]

where $F$ is a large constant and $\sum_{L}^{\text{odd}}$ only sums over every other layer (in a checkerboard pattern). If we are interested in the 1+1D duality (Sec. III.1), then we apply Eq. (18) and find that $H_{\text{inter}}^{\text{dual}}$ is dual to

[TABLE]

where $\sum_{\text{x-axis }\ell}^{y+z\text{ odd}}$ only sums over x-axis links $\ell$ with odd $y+z$. To describe the intermediate phase between X-cube order and YZ-plane 2+1D $Z_{2}$ QSL (intermediate #4 in Tab. 1), we can add $H_{\text{inter}}^{\text{1+1D}}$ to $\hat{H}_{\text{X-cube}}$ (Eq. (21)) and apply degenerate perturbation theory (DPT) using $H_{\text{inter}}^{\text{1+1D}}$ as the unperturbed Hamiltonian. This will generate a Hamiltonian which describes the intermediate phase. This process can be repeated to generate Hamiltonians for the possible intermediate phases for all four phase transitions shown in Fig. 1. We do not write down these Hamiltonians here, nor have we analyzed what kind of topological order is present in these models. However, we have calculated their ground state degeneracies, which we show in Tab. 1. This was done via a computer program Slagle (2017b) on finite lattices; we then extrapolated to arbitrary lattice sizes.