Twisted foliated fracton phases

TL;DR

This paper classifies various three-dimensional fracton models based on their foliated fracton order, revealing new twisted phases distinguished by their fractional excitations and subsystem symmetries.

Contribution

It identifies and constructs new twisted foliated fracton phases, expanding the understanding of fracton order beyond the well-known X-cube model.

Findings

Many known type I fracton models share the same X-cube foliated order.

New twisted fracton models are constructed with nontrivial order.

Twisted models exhibit distinct fractional excitations and subsystem symmetry properties.

Abstract

In the study of three-dimensional gapped models, two-dimensional gapped states should be considered as a free resource. This is the basic idea underlying the notion of `foliated fracton order' proposed in Phys. Rev. X 8, 031051 (2018). We have found that many of the known type I fracton models, although they appear very different, have the same foliated fracton order, known as `X-cube' order. In this paper, we identify three-dimensional fracton models with different kinds of foliated fracton order. Whereas the X-cube order corresponds to the gauge theory of a simple paramagnet with subsystem planar symmetry, the novel orders correspond to twisted versions of the gauge theory for which the system prior to gauging has nontrivial order protected by the planar subsystem symmetry. We present constructions of the twisted models and demonstrate that they possess nontrivial order by studying their fractional excitation contents.

Equations

$$e_A^2=e_B^2=1m_A^2=e_B{m}_B^2=e_A.$$

$$H_{TC}=-\sum _v\left(A_v+A_v^2+A_v^3\right)-\sum _p\left(B_p+B_p^2+B_p^3\right)$$

$$H_{2D}=-\sum _v\left(Q_v^A+Q_v^B\right)-\sum _p\left(O_p^A+O_p^B+\text{h.c.}\right),$$

$$O_p^A$$

$$O_p^B$$

$$\begin{array}{rl}{H}_{3D}=& -\sum _v\sum _{\sigma =x,y,z}\left(Q_{v,\sigma }^A+Q_{v,\sigma }^B\right)\\ & -\sum _c\left(O_c^A+O_c^B+\text{h.c.}\right),\end{array}$$

$${\left(f_{ijk}^A\right)}^2={\left(f_{ijk}^B\right)}^2={\left(e_{\mu \nu ,i}^A\right)}^2={\left(e_{\mu \nu ,i}^B\right)}^2=1.$$

$$\begin{array}{r}{\left(l_{\mu ,ij}^A\right)}^2=e_{\mu \nu ,i}^B\times e_{\rho \mu ,j}^B\\ {\left(l_{\mu ,ij}^B\right)}^2=e_{\mu \nu ,i}^A\times e_{\rho \mu ,j}^A\\ {\left(l_{\mu ,ij}^A\right)}^4={\left(l_{\mu ,ij}^B\right)}^4=1.\end{array}$$

$$l_x^A\times l_y^A\times l_z^A=l_x^B\times l_y^B\times l_z^B=1.$$

$$\begin{array}{rl}{m}_{\mu \nu ,i,i+1}^A& =l_{\mu ,ij}^A\times {\bar{l}}_{\mu ,i+1,j}^A=l_{\nu ,ki}^A\times {\bar{l}}_{\nu ,k,i+1}^A\\ m_{\mu \nu ,i,i+1}^B& =l_{\mu ,ij}^B\times {\bar{l}}_{\mu ,i+1,j}^B=l_{\nu ,ki}^B\times {\bar{l}}_{\nu ,k,i+1}^B,\end{array}$$

$$\begin{array}{rl}{\left(m_{\mu \nu ,i,i+1}^A\right)}^2& =e_{\mu \nu ,i}^B\times e_{\mu \nu ,i+1}^B\\ {\left(m_{\mu \nu ,i,i+1}^B\right)}^2& =e_{\mu \nu ,i}^A\times e_{\mu \nu ,i+1}^A.\end{array}$$

$$l_x^A\times l_y^A\times l_z^A=l_x^B\times l_y^B\times l_z^B=1.$$

$$\text{GSD}=2^{4L_x+4L_y+4L_z-6}.$$

$$I(A;B|C)=I(A;B;C;D|E)=\log (4).$$

$$\prod _i{e}_i=1,\prod _i{m}_{i,i+1}=1$$

$${\left(m_{i,i+1}^{A,B}\right)}^2=e_i^{B,A}\times e_{i+1}^{B,A}$$

$$\prod _i{m}_{i,i+1}^{B,A}=\prod _i{e}_i^{A,B}=1$$

$$m_i^2=e_{i-1}\times e_{i+1}.$$

$$\mathcal{L}=\frac{1}{4\pi }{K}_{IJ}{ϵ}^{\mu \nu \rho }{a}_{\mu }^I{\partial }_{\nu }{a}_{\rho }^J,$$

$${\theta }_l=\pi l^T{K}^{-1}l,$$

$${\theta }_{l{l}^{\prime }}=2\pi l^T{K}^{-1}{l}^{\prime }.$$

$$K=\bordermatrix{&&e_{1}&m_{1}&e_{2}&m_{2}&e_{3}&m_{3}&e_{4}&\cr&\ddots&&&&&&&&\cr&&0&2&\raisebox{0.75pt}{-}1&&&&&\cr&&2&0&&&&&&\cr&&\raisebox{0.75pt}{-}1&&0&2&\raisebox{0.75pt}{-}1&&&\cr&&&&2&0&&&&\cr&&&&\raisebox{0.75pt}{-}1&&0&2&\raisebox{0.75pt}{-}1&\cr&&&&&&2&0&&\cr&&&&&&\raisebox{0.75pt}{-}1&&0&\cr&&&&&&&&&\ddots}.$$

$$K^{-1}=\frac{1}{4}\bordermatrix{&&m_{0}&e_{1}&m_{1}&e_{2}&m_{2}&e_{3}&m_{3}&\cr&\ddots&&&&&&&&\cr&&0&&1&&&&&\cr&&&0&2&&&&&\cr&&1&2&0&&1&&&\cr&&&&&0&2&&&\cr&&&&1&2&0&&1&\cr&&&&&&&0&2&\cr&&&&&&1&2&0&\cr&&&&&&&&&\ddots}.$$

$$W=\bordermatrix{&&\tilde{e}_{1}&\tilde{m}_{1}&\tilde{e}^{A}&\tilde{m}^{A}&\tilde{e}^{B}&\tilde{m}^{B}&\tilde{e}_{2}&\tilde{m}_{2}&\tilde{e}_{3}&\cr&\ddots&&&&&&&&&&\cr e_{1}&&1&&&&&&&&&\cr m_{1}&&&1&&&&1&&&&\cr e_{2}&&&&1&&&&\raisebox{0.75pt}{-}1&&&\cr m_{2}&&&&&1&&&&&&\cr e_{3}&&\raisebox{0.75pt}{-}1&&&&1&&&&&\cr m_{3}&&&&&&&1&&&&\cr e_{4}&&&&&&&&1&&&\cr m_{4}&&\raisebox{0.75pt}{-}1&&&1&&&&1&&\cr e_{5}&&&&&&&&&&1&\cr&&&&&&&&&&&\ddots}.$$

$$W^{T}KW=\bordermatrix{&&\tilde{e}_{1}&\tilde{m}_{1}&\tilde{e}^{A}&\tilde{m}^{A}&\tilde{e}^{B}&\tilde{m}^{B}&\tilde{e}_{2}&\tilde{m}_{2}&\tilde{e}_{3}&\cr&\ddots&&&&&&&&&&\cr&&0&2&&&&&\raisebox{0.75pt}{-}1&&&\cr&&2&0&&&&&&&&\cr&&&&0&2&\raisebox{0.75pt}{-}1&0&&&&\cr&&&&2&0&0&0&&&&\cr&&&&\raisebox{0.75pt}{-}1&0&0&2&&&&\cr&&&&0&0&2&0&&&&\cr&&\raisebox{0.75pt}{-}1&&&&&&0&2&\raisebox{0.75pt}{-}1&\cr&&&&&&&&2&0&&\cr&&&&&&&&\raisebox{0.75pt}{-}1&&0&\cr&&&&&&&&&&&\ddots}.$$

$$K_{2D}^{-1}=\frac{1}{4}\bordermatrix{&\tilde{e}^{A}&\tilde{m}^{A}&\tilde{e}^{B}&\tilde{m}^{B}\cr&0&2&0&0\cr&2&0&0&1\cr&0&0&0&2\cr&0&1&2&0}.$$

$$m_i^2=e_{i-1}\times e_{i+1}$$

$$H_{\text{SSPT}}=-\sum _c\left({\tilde{O}}_c^A+{\tilde{O}}_c^B+\text{h.c.}\right)$$

$$U^{\alpha}_{P_{\mu\nu}}=\prod_{c\in P_{\mu\nu}}\sigma^{\alpha}_{x,c}\;\;\text{ with }\begin{array}[]{ll}\;\,\alpha=A,B\\ \mu\nu=xy,yz,zx\end{array}$$

$$H=-\sum _v\left(A_v^x+A_v^y+A_v^z\right)-\sum _c{B}_c$$

$${\log }_2\text{GSD}=2L_x+2L_y+2L_z-3.$$

Full text

I Introduction

The discovery of various “fracton” models Nandkishore and Hermele (2018); Vijay et al. (2015, 2016); Vijay (2017); Chamon (2005); Haah (2011); Yoshida (2013); Slagle and Kim (2018); Prem et al. (2018); Vijay and Fu (2017); Ma et al. (2017); Bravyi et al. (2011); Hsieh and Halász (2017); Halász et al. (2017); Slagle and Kim (2017a); Bulmash and Barkeshli (2018, 2019); Prem and Williamson (2019); Tian et al. (2018); Tian and Wang (2019); You et al. (2018) has greatly expanded our understanding of gapped phases in three-dimensional systems. A salient feature characterizing this set of models is the existence of gapped fractional point excitations with restricted mobility. Gapped fracton models111There are also gapless $U(1)$ fracton models which will not be addressed in this work (see Refs. Pretko, 2017a, b; Rasmussen et al., 2016; Pretko and Radzihovsky, 2018; Pretko, 2017c; Bulmash and Barkeshli, 2018; Ma et al., 2018; Gromov, 2018; Kumar and Potter, 2018; Pretko, 2018; Yan et al., 2019; Williamson et al., 2018b; Slagle et al., 2018). are divided into two major classes according to how the motion of point excitations is constrained: type I and type II. In type II models, fractional point excitations can only move in coordination as a set and individually they cannot move at all. These excitations are said to be immobile and are called ‘fractons’. In type I models, on the other hand, apart from fracton excitations, there can also exist lineons and planons – fractional excitations which can move by themselves within a plane or along a line. The restricted mobility of the point excitations leads to various new features in the fracton models: a slow thermalization process,Bravyi and Haah (2011); Prem et al. (2017); Pai et al. (2018); Pai and Pretko (2019); Sala et al. (2019) stable extensive ground state degeneracy, unusual entanglement scaling, Shirley et al. (2018a); Shi and Lu (2018); Ma et al. (2018); He et al. (2018) etc.Yan (2018); Bulmash and Iadecola (2018); Pai and Hermele (2019); Brown and Williamson (2019); Williamson et al. (2018a); You (2019); Kubica and Yoshida (2018)

Among the type I fracton models, we have found that many of them have a hidden ‘foliation’ structure and are said to have ‘foliated fracton order’ (FFO).Shirley et al. (2017, 2018a) That is, starting from a model with a larger system size, we can apply a finite depth local unitary transformation and map the model to a smaller system size together with decoupled layers of 2D gapped states, as illustrated in Fig. 1(a). As there should be no fundamental change in the order of the system simply due to the change in system size, we should think of the 2D gapped states as free resources in the study of these 3D fracton models even though the 2D gapped states can have highly nontrivial topological order of their own. Correspondingly, we define two foliated fracton models to have the same ‘foliated fracton order’ if they can be related through a finite depth local unitary transformation upon the addition of decoupled stacks of 2D layers of gapped states, as shown in Fig. 1(b).

Using this definition, we have found that many of the type I fracton models with a foliation structure actually have the same foliated fracton order. In particular, we have shown explicitly that the spin checkerboard model, the Majorana checkerboard model, and the semionic X-cube model all have the same FFO as the X-cube model (or multiple copies of it). Shirley et al. (2019a); Wang et al. (2019); Shirley et al. (2018b) The untwisted string-membrane-net model discussed in Ref. Slagle et al., 2019 was also shown to be equivalent to the X-cube model. As the X-cube modelVijay et al. (2016) can be obtained by gauging the intersecting planar subsystem symmetries of a trivial 3D paramagnet,Vijay et al. (2016); Shirley et al. (2019b) the X-cube FFO is considered to be untwisted. It is similar to the toric code model as an untwisted $Z_{2}$ gauge theory which can be obtained by gauging the global $Z_{2}$ symmetry of a trivial 2D paramagnet. It is known that 2D $Z_{2}$ gauge theory can also be ‘twisted’ where the gauge flux becomes a semionic excitation. It can be obtained from gauging the 2D symmetry protected topological order with $Z_{2}$ symmetry as shown in Ref. Levin and Gu, 2012. It is then natural to ask whether there exists twisted FFO.

In this paper, we identify three-dimensional fracton models with a ‘twisted’ FFO. That is, these models have an FFO that is different from that of the X-cube model. Moreover, they can be obtained by gauging a 3D model with subsystem planar symmetries that is not a trivial paramagnet. In other words, the ungauged model has (strong) symmetry protected topological order with subsystem planar symmetries. Note that although twisted fracton models have already appeared in the literature, You et al. (2018a); Song et al. (2018); Devakul et al. (2018) they have not been studied in terms of their foliated fracton order. We discuss two (sets of) examples in detail. One is 3-foliated, meaning that we can decouple 2D topological layers in three different directions using finite depth local unitary transformations. The X-cube model is also 3-foliated in this sense and we can consider this new model as the twisted version of two copies of the X-cube model. The other example is 1-foliated, meaning that we can only decouple 2D topological layers in one direction from the model. The untwisted version of a 1-foliated model with $Z_{2}$ symmetries would simply be a decoupled stack of 2D toric codes.

The paper is organized as follows: In section II, we discuss the 3-foliated model by presenting the construction of the model, demonstrating its foliation structure and then showing that its FFO is different from that of the X-cube model. In section III, we do the same for the 1-foliated model. We discuss in section IV how to ‘ungauge’ the models into models with subsystem symmetry protected topological order before summarizing in section V.

II Twisted 3-Foliated Model

In this section, we describe a model that is foliated in the $x$, $y$, and $z$ directions by layers of a twisted 2D $Z_{2}\times Z_{2}$ gauge theory. We will see that its foliated fracton order (FFO) is twisted in the sense that its FFO is distinct from that of the X-cube model or copies of it. (A brief review of the X-cube model is given in Appendix A.) Ungauging this fracton model results in a paramagnetic model with (strong) subsystem symmetry protected topological (SSPT) order under 3 sets of intersecting planar subsystem symmetries.

The model is constructed by strongly coupling intersecting layers of a set of 3 perpendicular stacks of twisted 2D $Z_{2}\times Z_{2}$ gauge theories, in a manner akin to the construction of the X-cube and semionic X-cube models from stacks of 2D toric code layers and 2D double semion layers respectively. These constructions are discussed in Refs. Ma et al., 2017 and Vijay, 2017. Like the semionic X-cube model, the 3-foliated model constructed in this section belongs to the class of exactly solvable twisted fracton models considered in Ref. Song et al., 2018. Here, we are able to extend our understanding by studying the model through the lens of the coupled layer construction and as an FFO. Unlike the semionic X-cube model, this $Z_{2}\times Z_{2}$ model has twisted FFO; thus, there is a distinction between a fracton model being twisted in the sense of Ref. Song et al., 2018, and a model having twisted FFO.

II.1 Model Construction

II.1.1 2D $Z_{2}\times Z_{2}$ twisted gauge theory

First, we briefly review the properties of $Z_{2}\times Z_{2}$ twisted gauge theories in 2D, and describe an exactly solvable model for one such theory. Twisted gauge theories may be thought of as Hamiltonian realizations of 2+1d Dijkgraaf-Witten models,Dijkgraaf and Witten (1990) or as the result of gauging global symmetries in paramagnets with non-trivial symmetry-protected topological (SPT) order.Levin and Gu (2012) For $Z_{2}\times Z_{2}$ symmetry, there are $2^{3}=8$ distinct SPT phases in 2D, corresponding to the 8 elements of $H^{3}(Z_{2}\times Z_{2},U(1))$.Chen et al. (2013) They are characterized by the topological invariants $N_{1}$, $N_{2}$, and $N_{12}$, each of which takes values 0 or 1. Upon gauging, the exchange statistics of the gauge fluxes are given by $i^{N_{1}}$ and $i^{N_{2}}$, whereas the braiding statistics between the two fluxes is $i^{N_{12}}$. In all cases the statistics between gauge charge and corresponding gauge flux is $-1$.Wang and Levin (2015)

Here, we will focus on the twisted gauge theory obtained from the SPT phase with $N_{12}=1$ and $N_{1}=N_{2}=0$. In this case, the elementary gauge charges $e_{A}$ and $e_{B}$ and bosonic gauge fluxes $m_{A}$ and $m_{B}$ obey the following fusion rules:

[TABLE]

Thus, as an intrinsic topological order, this theory is equivalent to the $Z_{4}$ toric code, with $m_{A}$ and $m_{B}$ mapping onto the $\tilde{e}$ and $\tilde{m}$ $Z_{4}$ anyons of the $Z_{4}$ toric code, respectively.

There is a convenient isomorphism between the $Z_{4}$ clock and shift algebra and the two qubit operator algebra, $\tilde{Z}\to X^{A}S^{B}$, $\tilde{X}\to X^{B}CZ^{AB}$, where $\tilde{Z}$ and $\tilde{X}$ are the clock and shift generators of the $Z_{4}$ operator algebra with $\tilde{Z}\tilde{X}=i\tilde{X}\tilde{Z}$, $Z$ and $X$ are the clock and shift (Pauli) operators of the $Z_{2}$ algebra, $S$ is the one-qubit phase gate $\text{diag}(1,i)$, $CZ$ is the two-qubit controlled-Z operator $\text{diag}(1,1,1,-1)$, and $A$ and $B$ label the two qubits. Applied to the $Z_{4}$ toric code degrees of freedom, this mapping naturally allows one to write the $Z_{2}\times Z_{2}$ twisted gauge theory as a $Z_{2}\times Z_{2}$ string-net model, such that the gauge charges correspond to violations of the plaquette terms.

In particular, consider the $Z_{4}$ toric code Hamiltonian on a square lattice:

[TABLE]

where $A_{v}=\tilde{Z}^{\dagger}_{1}\tilde{Z}^{\dagger}_{5}\tilde{Z}_{6}\tilde{Z}_{2}$ and $B_{p}=\tilde{X}_{1}\tilde{X}_{2}\tilde{X}^{\dagger}_{3}\tilde{X}^{\dagger}_{4}$ per Fig. 2. The $A_{v}^{2}$ and $B_{p}^{2}$ terms are redundant, but we keep them in the Hamiltonian so that the transformed Hamiltonian has a convenient correspondence with the string-net formulation.Levin and Wen (2005)

After mapping to qubit degrees of freedom and shifting qubit $A$ downward and to the left by half a lattice spacing [as shown in Fig. 2(a)], $H_{TC}$ is transformed into the $Z_{2}\times Z_{2}$ twisted gauge theory Hamiltonian

[TABLE]

where $Q^{\mu}_{v}=\prod_{l\in v}Z^{\mu}_{l}$, and (see Fig. 2)

[TABLE]

In particular, $A_{v}\to O_{p}^{A}$, $A_{v}^{3}\to\left(O_{p}^{A}\right)^{\dagger}$, $B_{p}\to O_{p}^{B}$, $B_{p}^{3}\to\left(O_{p}^{B}\right)^{\dagger}$, $A_{v}^{2}\to Q_{v}^{A}$, and $B_{p}^{2}\to Q_{v}^{B}$. Note that $\left(X^{B}CZ^{AB}\right)^{\dagger}=X^{B}CZ^{AB}Z^{A}$.

As this transformation is an exact mapping, it is obviously possible to carry through the following construction in terms of the original $Z_{4}$ degrees of freedom. As we will see however, the $Z_{2}\times Z_{2}$ degrees of freedom provide a more natural language to analyze the emergent fracton order.

II.1.2 Coupled layers construction

The construction of the 3-foliated fracton model is a straightforward generalization of the construction of the X-cube and semionic X-cube models in Refs. Ma et al., 2017. We first start with 3 mutually perpendicular intersecting stacks of the $Z_{2}\times Z_{2}$ twisted gauge theory model $H_{2D}$, oriented as in Fig. 3. Recall that $H_{2D}$ contains 2 qubit degrees of freedom ($A$ and $B$) on each edge of a square lattice. Each edge of the 2D layers coincides with another edge from an orthogonal layer to form a cubic lattice, with each edge containing 4 qubits. Then, couplings of the form $Z^{A}Z^{A}$ and $Z^{B}Z^{B}$ between qubits on the same edge are added to the Hamiltonian.

In the strong coupling limit, the four qubits at each edge merge into two. The following effective Hamiltonian emerges at lowest order in perturbation theory:

[TABLE]

where $v$ runs over vertices of the cubic lattice and $c$ runs over the elementary cubes. $Q^{\mu}_{v,\sigma}$ are vertex terms equal to products of Pauli $Z^{\mu}$ operators over the links adjacent to $v$ in the plane normal to $\sigma$. The cube operators $O^{A}_{c}$ and $O^{B}_{c}$ are depicted graphically in Fig. 3. The terms of $H_{3D}$ are mutually commuting and unfrustrated and thus the model is exactly solvable. It bears striking similarity to (two copies of) the X-cube model: the vertex terms are identical, and the cube terms are similar in that they involve products of Pauli $X$ operators over the edges of the cube. However, they contain additional phase factors not present in the X-cube terms.

As in the X-cube model, excitations of the vertex constraints are lineons whereas excitations of the cube terms are fractons. Lineons are created at the endpoints of open rigid string operators, whereas fractons are created at the corners of membrane operators. Examples of these operators are given in the discussion of interferometric operators in Sec. II.3.2. Like the X-cube model, planons also exist as fracton dipoles and lineon dipoles, as will be discussed in detail below.

II.2 Fractional excitations

In the intermediate coupling regime, the transition to the strong-coupling phase can be thought of as a condensation of $A$ and $B$ type charge loops; correspondingly the ground state of $H_{3D}$ may be viewed as a condensate of charge loops. This mechanism has been studied in detail and dubbed $p$-string condensation in Ref. Ma et al., 2017. The structure of excitations in the condensed phase can be understood in terms of the degrees of freedom of the pre-condensed stacks of twisted $Z_{2}\times Z_{2}$ gauge theories. Similar to the case of the X-cube model discussed in Refs. Ma et al., 2017, the 2D gauge charges of the original decoupled stacks fractionalize into fracton dipoles (a pair of adjacent fractons whose axis is normal to the 2D layer), and remain as $Z_{2}$ planons. These planons will be labelled $e^{A}_{\mu\nu,i}$ and $e^{B}_{\mu\nu,i}$ where $\mu\nu,i$ refers to the plane of mobility ($\mu$ and $\nu$ the planar axes and $i$ the coordinate in the normal direction). In the charge loop picture, individual fractons correspond to endpoints of open charge strings above the condensate. They will be denoted as $f^{A}_{ijk}$ and $f^{B}_{ijk}$, where $ijk$ denotes spatial location, and likewise inherit $Z_{2}$ fusion rules:

[TABLE]

As in the X-cube coupled layers construction, individual gauge fluxes of the original stacks are confined upon condensation due to their statistical interaction with the charge loops. However, composites of an $A$ ($B$) flux and an $A$ ($B$) anti-flux in orthogonal planes have trivial statistics with the charge loops, and thus survive the condensation. These composites become $A$ and $B$ type lineons of the condensed phase, labelled as $l^{A}_{\mu,ij}$ and $l^{B}_{\mu,ij}$ with $\mu$ the axis of mobility and $i$ and $j$ the normal coordinates. By convention $l^{A}_{\mu,ij}$ ($l^{B}_{\mu,ij}$) consists of a flux in the $\mu\nu$ plane and an anti-flux in the $\rho\mu$ plane. They inherit the fusion rules from the 2D gauge fluxes, and therefore obey:

[TABLE]

In these equations, the fracton dipoles’ planes of mobility intersect along the lineon axis. There are also triple fusion rules between intersecting lineons along orthogonal axes (coordinate labels have been suppressed):

[TABLE]

Whereas individual lineons are restricted to move along a line, adjacent lineon anti-lineon pairs, called lineon dipoles, are free to move in a plane normal to the axis of separation, and are hence planons. This is because lineons arise as bound states of flux anti-flux pairs in orthogonal planes. A lineon dipole therefore contains four original flux (or anti-flux) excitations. However, the flux anti-flux pair in the plane shared by the two lineons annihilate one another, leaving behind a flux anti-flux pair in adjacent parallel planes. Lineon dipoles will be denoted $m^{A}_{\mu\nu,i,i+1}$ and $m^{B}_{\mu\nu,i,i+1}$ where $\mu\nu$ refers to the plane of mobility and $i$ and $i+1$ are the coordinates in the normal direction of the parallel planes containing the flux and anti-flux respectively. The following fusion rules hold by definition:

[TABLE]

where $\bar{l}$ refers to the anti-lineon of $l$. Combining Eq. (7) and Eq. (9) yields the rules

[TABLE]

The statistics of excitations in the condensed phase can also be inferred from the anyon statistics of the decoupled stacks. In particular, the fracton dipole $e^{A}_{\mu\nu,i}$ ($e^{B}_{\mu\nu,i}$) exhibits a $-1$ braiding statistic when wound around type $A$ ($B$) lineons mobile within the dipole’s plane of movement. In particular, these lineons are $l^{A}_{\nu,ij}$ ($l^{B}_{\nu,ij}$) and $l^{A}_{\mu,ji}$ ($l^{B}_{\mu,ji}$). Moreover, coplanar lineons of opposite species $l^{A}_{\mu,ij}$ and $l^{B}_{\nu,ki}$ inherit the $i$ braiding statistic between gauge fluxes $m_{A}$ and $m_{B}$; thus they exhibit an $i$ statistical phase upon crossing. This property, along with the lineon fusion rules, are the essential features that distinguish the twisted 3-foliated model from the untwisted version, i.e. two copies of the X-cube model.

II.3 Foliation structure

In this section, we first show that the model described in the last section indeed has a foliated fracton order. That is, one can decouple 2D topological layers out of the model while shrinking the system size as shown in Fig. 1 (a). Then we are going to look at some of the universal quantities of foliated fracton orders, including the quotient super-selection sectors and the entanglement signatures that we discussed in Refs. Shirley et al., 2018b and Shirley et al., 2018a. It turns out that this model is trivial (the same as two copies of the X-cube model) in both aspects. However, it is not equivalent to two copies of the X-cube model as an FFO, which we will show in Sec. II.4.

II.3.1 Resource layers

In this section, we demonstrate the 3-foliated structure of the model. We show that resource layers consisting of bilayer 2D $Z_{2}\times Z_{2}$ twisted gauge theories can be decoupled from the model in all three directions. Rather than finding an exact local unitary transformation, we arrive at this conclusion by examining the structure of fractional excitations in an $L_{x}\times L_{y}\times L_{z}$ size 3D model, and find that it can be decomposed into two parts: one corresponding to a reduced $L_{x}\times L_{y}\times(L_{z}-2)$ size 3D model, and the other corresponding to two layers of the twisted gauge theory described by $H_{2D}$. That is, the superselection sectors of the larger 3D model are identical to those of the smaller 3D model together with the decoupled 2D layers. We may then conclude the presence of such a foliation structure.

In gapped abelian phases, the superselection sectors form an abelian group under fusion. Decomposing this structure therefore amounts to finding a generating set of the fusion group which can be bipartitioned into sets $A$ and $B$ such that there are no statistical interactions between sectors of $A$ and sectors of $B$.

For the model in question, $S$ contains fractons, lineons, and planons. However, the elementary planons are either fracton dipoles or lineon dipoles (lineon anti-lineon pairs). Therefore, fusion with the appropriate planon effectively transports lineons or fractons in their directions of immobility. Hence, a generating set of $S$ need only one lineon of each type in each direction, one fracton of each type, and a generating set of the planon subgroup $P\leq S$ (i.e. the subgroup of $S$ generated by the set of all planons), which decomposes as $P=P_{xy}\times P_{yz}\times P_{zx}$ for the three different planes of mobility. This phenomenon also occurs in all of the stabilizer code models with FFO that have been previously studied. In fact, this observation is the basis of the notion of quotient superselection sectors (QSS), which are elements of the quotient group $Q=S/P$ to be discussed below.

Suppose we wish to disentangle a resource layer in the $z$ direction from the twisted 3-foliated model. Due to the above observation, a decomposition $S=S_{2D}\times S^{\prime}$, where $S_{2D}$ represents a single 2D resource layer and $S^{\prime}$ is the reduced 3D model, amounts to a decomposition $P=P_{2D}\times P^{\prime}=\left(P_{2D}\times P_{xy}^{\prime}\right)\times P_{yz}\times P_{zx}$, such that $P_{2D}$ has no statistical interaction with $P^{\prime}$. Moreover, $P_{2D}$ must have trivial interactions with the generating lineons and fractons. However, these generators can always be chosen to lie away from the support of the $P_{2D}$ string operators; thus, this latter condition is essentially vacuous.

Let us now consider $P_{xy}$, the subgroup of $S$ consisting of planons mobile in the $x$ and $y$ directions. A generating set of $P_{xy}$ is given by the set of elementary (minimally separated) fracton and lineon dipoles with $z$-oriented dipolar axis. It is possible to find an equivalent generating set that decouples into two subsets: one generates $P_{xy}^{\prime}$, a reduced version of $P_{xy}$; the other generates $P_{2D}$, which corresponds to two copies of the 2D $Z_{2}\times Z_{2}$ twisted gauge theory modeled by $H_{2D}$. To illustrate this decomposition, it is convenient to use a graphical notation, as shown in Fig. 4. In Fig. 4(a), (part of) the generating set of elementary dipoles is depicted. Fig. 4(b) contains an equivalent but different generating set. In this set, the quasiparticles represented by rows 5-12 are completely decoupled from the remaining planons, in the sense that they form a closed group under fusion and have trivial braiding statistics with the other planons. These quasiparticles represent a generating set of the anyon sectors of two copies of the $Z_{2}\times Z_{2}$ twisted gauge theory, i.e. a bilayer (rows 5-8 and rows 9-12). The remaining planons constitute a reduced version of the original planon group with two fewer lattice spacings in the $z$ direction.

Importantly, this mapping of generating planons preserves the locality of the excitations in the $z$ direction. In other words, each element of the generating set moves within a finite region in $z$ before and after the mapping. Therefore, we expect that this mapping of excitations can be realized by a finite depth local unitary transformation with support in the vicinity of the decoupled resource bilayer.

Having established the foliation structure in the 3-foliated model, we now ask if it has the same FFO as the X-cube model (or copies of it). As defined in Refs. Shirley et al., 2017 and Shirley et al., 2018a, two gapped models have the same foliated fracton order (FFO) if they can be related by a local unitary transformation upon the possible addition of 2D topological order resource states. While this is a rather coarse equivalence relation, previous works have identified the structure of QSS and interferometric statistics, as well as the entanglement signatures discussed prior, as universal characteristics of FFO.Shirley et al. (2018a, b) As we are going to see in section II.3.2 and II.3.3, based on these properties alone it is plausible that the 3-foliated model has the same FFO as two copies of the X-cube model. However, as we are going to show in section II.4, the 3-foliated model actually has a different FFO from two copies of the X-cube model. The QSS and entanglement signature hence provide an insufficient characterization of the universal properties of a foliated fracton phase.

II.3.2 Quotient superselection sectors and interferometric statistics

Consider the QSS fusion group $Q=S/P$. To reiterate, the essential idea behind QSS is that by modding out the planon subgroup $P$, we obtain a finite group which is characteristic of the foliated fracton order of a given model. Since lineon and fracton dipole sectors belong to $P$ for the twisted 3-foliated model, it follows that all lineon superselection sectors $l^{A}_{\mu,ij}$ ($l^{B}_{\mu,ij}$) belong to one quotient sector, denoted $l^{A}_{\mu}$ ($l^{B}_{\mu}$). Moreover all fracton sectors $f^{A}_{ijk}$ ($f^{B}_{ijk}$) belong to a single quotient sector, denoted $f^{A}$ ($f^{B}$). These quotient sectors generate the entire group $Q$.

However, lineon and fracton quotient sectors also obey some relations. First, since $e^{A}_{\mu\nu,i}\times e^{A}_{\rho\mu,j}$ and $e^{B}_{\mu\nu,i}\times e^{B}_{\rho\mu,j}$ belong to $P$, the lineon fusion rules (7) imply that $\left(l^{A}_{\mu}\right)^{2}=\left(l^{B}_{\mu}\right)^{2}=1$ as quotient sectors. In other words, the lineon quotient sectors obey $Z_{2}$ fusion rules while the lineon superselection sectors obey $Z_{4}$ fusion rules. Second, the lineon triple fusion rules are inherited by the quotient group as

[TABLE]

Finally, the fractons sectors obey $\left(f^{A}\right)^{2}=\left(f^{B}\right)^{2}=1$. Therefore, altogether $Q\cong\left(Z_{2}\right)^{\times 6}$, with the generators $f^{A}$, $f^{B}$, $l^{A}_{x}$, $l^{B}_{x}$, $l^{A}_{y}$, and $l^{B}_{y}$. This QSS structure is isomorphic (in terms of fusion and particle mobility) to that of two copies of the X-cube model, one corresponding to each of the $A$ and $B$ sectors of $Q$. Recall that the X-cube model has QSS group $\left(Z_{2}\right)^{\times 3}$ with generators $f$, $l_{x}$, and $l_{y}$, and triple fusion rule $l_{x}\times l_{y}\times l_{z}=1$.

Interferometric operators for foliated orders, as introduced in Ref. Shirley et al., 2018b, are unitary operators with support outside the region $R$, where a point excitation is located, that yield nontrivial statistical phases when acting on excitations belonging to nontrivial elements of $Q$, but act as the identity on excitations in $P$. As discussed in Ref. Shirley et al., 2018b, for the X-cube model, there are 8 classes of such operators, which have a $Z_{2}\times Z_{2}\times Z_{2}$ group structure. They include a wireframe operator $W$ which yields a $-1$ phase on the quotient sector $f$, and cylindrical membrane operators $M_{x}$, $M_{y}$, and $M_{z}$. The operator $M_{x}$ yields a $-1$ phase on the $l_{y}$ and $l_{z}$ sectors, and similarily for $M_{y}$ and $M_{z}$.

In the twisted 3-foliated model, the structure of interferometric operators is identical to that of two copies of the X-cube model, in terms of the geometry of the operators and their statistical interactions with the QSS. In particular, there are operators $W^{A}$, $W^{B}$, $M_{x}^{A}$, $M_{x}^{B}$, $M_{y}^{A}$, $M_{y}^{B}$, $M_{z}^{A}$, and $M_{z}^{B}$. The microscopic form of these operators may be computed by taking products of all the Hamiltonian terms of one kind within a large cubic region: the wireframe operators $W^{A}$ and $W^{B}$ correspond to products of cube operators $O^{A}_{c}$ and $O^{B}_{c}$, whereas the membrane operators correspond to products of the vertex terms. Thus, the membrane operators are simply products of Pauli $Z^{A}$ or $Z^{B}$ operators over the support of the membrane, as in (two copies of) the X-cube model, whereas the wireframe operators are more complicated.

The rigid string and membrane operators, which create and transport lineons and fractons, have the identical form as these interferometric operators away from the excitations. The statistical interactions between interferometric operators and QSS can be verified by considering the commutation relations of these microscopic operators. One may also view the interferometric operators as planon loop operators for lineon or fracton dipoles with a macroscopic dipolar length.

II.3.3 Ground state degeneracy and entanglement signatures

To efficiently calculate the entanglement properties of the 3-foliated model, we consider a $Z_{4}$ CSS stabilizer code formulation of the model. That is, the Hamiltonian can be expressed as a sum of products of either $\tilde{Z}$ (the clock operator) or $\tilde{X}$ (the shift operator) where all terms in the Hamiltonian commute with each other and each term has eigenvalue $-1$ in the ground state. This form of Hamiltonian is useful for doing computations, and will allow us to efficiently calculate ground state degeneracy and entanglement entropy. In Appendix C, we will also express this model in the string-membrane-net and foliated field theory formulations.

To obtain a CSS version of the model, we can repeat the coupled layer construction from Sec. II.1.2, but continue using the $Z_{4}$ clock and shift operators instead of mapping to pairs of qubits. The coupled layer construction was performed by adding Pauli $Z^{A}Z^{A}$ and $Z^{B}Z^{B}$ terms to couple the $Z_{2}\times Z_{2}$ twisted gauge theory layers together. The $Z^{A}$ and $Z^{B}$ operators are written in terms of $\tilde{Z}^{2}$ and $\tilde{X}^{2}$, as in Fig. 2(a). Thus, the $Z^{A}Z^{A}$ term that couples $Z_{2}\times Z_{2}$ twisted gauge theory layers is mapped back to a $\tilde{Z}^{2}\tilde{Z}^{2}$ term to couple $Z_{4}$ toric code layers together. Unmapping the $Z^{B}Z^{B}$ term is similar, although note that the $\tilde{X}^{2}$ operator is not on the same edge as the $Z^{B}$ operator. Therefore, the $Z^{B}Z^{B}$ term is mapped back to a $\tilde{X}^{2}\tilde{X}^{2}$ operator, but where each $\tilde{X}^{2}$ is on a different link. The strong coupling limit is described by the CSS code Hamiltonian in Fig. 5.

Since the model is a stabilizer code, we can efficiently calculate its ground state degeneracy and entanglement entropy (see Appendix B for details). The ground state degeneracy of an $L_{x}\times L_{y}\times L_{z}$ system with periodic boundary conditions is

[TABLE]

Two-dimensional topological orders can be characterized by their topological entanglement entropy. Kitaev and Preskill (2006); Levin and Wen (2006) Ref. Shirley et al., 2018a discussed a generalization for foliated fracton orders given by the entanglement quantities $I(A;B|C)$ and $I(A;B;C;D|E)$ computed from subsystems with the wireframe geometries shown in Fig. 6. For the 3-foliated Hamiltonian (Eq. (5)), we find that

[TABLE]

These entanglement signatures, as well as the ground state degeneracy, are equivalent to that of two copies of the X-cube model.222In Ref. Shirley et al., 2018a, logarithms were evaluated in base 2. With this convention, the entanglement quantities in Eq. (12) are $I(A;B|C)=I(A;B;C;D|E)=\log_{2}(4)=2$. The X-cube model has $I(A;B|C)=I(A;B;C;D|E)=\log_{2}(2)=1$.

II.4 Twisted foliated fracton order

While the 3-foliated model appears the same as two copies of the X-cube model in terms of QSS and entanglement signatures, they actually have different FFO. In this section we will demonstrate this difference in two separate ways.

II.4.1 Lineon fusion rules

First, we will show that the $Z_{4}$ fusion rules of the lineon superselection sectors preclude a transformation to two copies of the X-cube model through local unitary and addition of 2D layers. It will be helpful to establish some terminology. A superselection sector that is a fusion product of planons in orthogonal planes, such that the mobility is restricted to the line of intersection of the two planes, will be referred to as a superficial lineon. Conversely, a lineon sector that cannot be decomposed as the fusion product of two planons, is referred to as an intrinsic lineon.Prem et al. (2018) While intersecting stacks of decoupled 2D topological orders exhibit superifical lineon superselection sectors, only truly fractonic models host intrinsic lineon excitations.

The key to the argument is that all of the intrinsic lineons in the twisted 3-foliated model are order 4 under fusion (although they square to superficial lineons hence the QSS has order 2), whereas the X-cube model contains intrinsic lineons of order 2. By adding stacks of 2D topological orders, it possible to modify the superselection sector group to include new intrinsic lineons of a higher order than the already existing intrinsic lineons. However, the fusion rules of the original intrinsic lineons are immutable, and moreover it is not possible to create a new intrinsic lineon of a lower order than the already existing sectors. Therefore, even after the free addition of 2D topological order resource states, the twisted 3-foliated model can never contain intrinsic lineons of order 2. Conversely, the X-cube model, and any number of copies of it, will always retain such a intrinsic lineon. Thus, the two models must have different FFO.

II.4.2 Redundancies among planons

Another way to see that the FFO of the 3-foliated model is different from that of two copies of the X-cube model is by looking at the planons. In fact, this can be a useful and generic way to study foliated fracton models. In the following, we are going to show that by examining the planons, we can deduce, first, that the X-cube model is different from a stack of 2D layers and secondly, that the 3-foliated model is different from the X-cube model (or 2 copies of it).

Consider a dimensional reduction procedure from a 3D model to a 2D model where the $x$ and $y$ directions remain infinite while the $z$ direction is made finite. Such a ‘compactification’ process has been used in Ref. Dua et al., 2019 to study fracton models. We consider the situation where the system has periodic boundary condition in all three directions. As the model is now finite in the $z$ direction, any string operator that extends around the $z$ direction becomes finite and can be added to the Hamiltonian. The ground state degeneracy is reduced and the model becomes a 2D model with anyons moving in the 2D plane. Here, we consider what happens upon this compactification process in three different fracton models: a decoupled stack of 2D layers, the X-cube model and the 3-foliated model.

We start with a decoupled stack of 2D layers in the $xy$ plane. In the 3D model, there is no string operator in the $z$ direction, therefore after dimensional reduction no extra term can be added. All the planons in the $xy$ planes survive the dimensional reduction. The number of planons grows exponentially with the height of the system in the $z$ direction. We can choose a generating set of all the planons by choosing a generating set for each plane. Such a generating set satisfies the following properties:

• •

Each element in the generating set is constrained to move within a finite segment in $z$ as they come from the 2D layers. We say that the generator planons are ‘local’.

• •

All other planons that are local can be generated by a subset of the generators that are within a finite distance in $z$. We say that the generating set is ‘locally complete’.

• •

Moreover, we can make sure the full generating set is not redundant. That is, no element in the generating set (or copies of it) can be generated by other elements in the generating set.

For the X-cube model and the 3-foliated model, these properties can no longer be satisfied at the same time.

Now we consider the X-cube model. A brief review of the X-cube model is given in Appendix A. Upon dimensional reduction in the $z$ direction, the string operators in the $z$ direction can be added to the Hamiltonian. Among all the fractional excitations, only the planons in the $xy$ planes survive the dimensional reduction procedure and we can choose a generating set for them consisting of the fracton dipoles $e_{i}$ centered around plane $i$ and the lineon dipoles $m_{i,i+1}$ living across planes $i$ and $i+1$. Such a generating set is local and locally complete as we defined above. However, it is redundant as the product of all fracton dipoles and the product of all lineon dipoles are both trivial anyons.

[TABLE]

That is, there exists global constraints among the planons. These global constraints cannot be removed without violating the ‘locally complete’ condition. If we remove $e_{1}$ and $m_{1}$ from the generating set, the set is no longer redundant, but $e_{1}$ and $m_{1}$ can not be locally generated. Therefore, the X-cube model is different from a stack of 2D layers.

Finally we turn to the 3-foliated model and see how it is different from both the stack of 2D layers and the X-cube model. Upon dimensional reduction, all other superselection sectors are removed except planons in the $xy$ plane, which are the fracton dipoles $e^{A,B}_{i}$ and lineon dipoles $m^{A,B}_{i,i+1}$. The $e^{A,B}_{i}$’s and $m^{A,B}_{i,i+1}$ sectors form a locally complete generating set, but it is highly redundant. First, there are local redundancies of the form

[TABLE]

Moreover, there are global redundancies of the form

[TABLE]

The global redundancies are similar to that of the X-cube, but the local ones show that the 3-foliated model is different from the X-cube. Note that it is possible to have local redundancy in a locally complete generating set of the X-cube model. For example, if besides all the $e_{i}$s and and $m_{i,i+1}$s we add $\psi_{i}=e_{i}\times m_{i,i+1}$ to the generating set, it will have a local redundancy. However, such local redundancies can be locally removed. That is, if we use the relation $\psi_{i}=e_{i}\times m_{i,i+1}$ and eliminate $\psi_{i}$ from the generating set, we can remove the redundancy. On the other hand, this is not true for the local redundancies in the 3-foliated model. In the 3-foliated model, we can start from the redundancy relation $\left(m^{A,B}_{1,2}\right)^{2}=e^{B,A}_{1}\times e^{B,A}_{2}$ and remove it by eliminating $e^{B,A}_{2}$ from the generating set. Next, we move on to eliminate $e^{B,A}_{3}$ from the generating set using the redundancy relation $\left(m^{A,B}_{1,2}\times m^{A,B}_{2,3}\right)^{2}=e^{B,A}_{1}\times e^{B,A}_{3}$. We can keep doing this, but the redundancy relation that we need to use involves more and more $m$ sectors, and eventually it becomes a non-local relation. We say that the local redundancy relations cannot be locally removed. In fact, a locally complete generating set always has to contain a finite density of $e$ particles and all the $m$ particles, therefore it is always redundant and the redundancy cannot be removed locally. Because of the existence of redundancy relations, especially local redundancy relations that cannot be locally removed, the 3-foliated model is different from both the stack of 2D layers and the X-cube model.

III Twisted 1-foliated model

In this section, we discuss a model which is non-trivially 1-foliated. That is, growing the model in the $z$ direction requires the addition of 2D topological order resource layers ($Z_{2}\times Z_{2}$ twisted gauge theories for the model we study), whereas growing the model in the $x$ or $y$ directions simply requires product state resources. At the same time, the model is not local unitarily equivalent to a decoupled stack of 2D topological orders. Nonetheless, all of the fractional excitations of the model are planons, which are mobile in the $xy$ directions; upon compactification in the $z$ direction,333I.e., a dimensional reduction from a 3D system to a 2D system with a large unit cell. the model reduces to a ‘giant’ 2D topological order where the number of superselection sectors grows exponentially with the original height in the $z$ direction.

III.1 Model construction

III.1.1 Boson condensation

The model is constructed by condensing bosons in a decoupled stack of 2D $Z_{2}\times Z_{2}$ twisted gauge theories (equivalently a stack of $Z_{4}$ toric codes, as discussed in Sec. II.1.1), stacked in the $z$ direction. The quasiparticle sectors of the stack consist of $Z_{2}\times Z_{2}$ gauge charges $e^{A}_{i}$ and $e^{B}_{i}$ and gauge fluxes $m^{A}_{i}$ and $m^{B}_{i}$. Composites of gauge charges in neighboring layers, $e^{A}_{i}e^{B}_{i-1}$, are then condensed to yield a new phase, whose fractional excitations can be understood in the conventional framework of 2D boson condensation in topological phases.Burnell (2018); Bais and Slingerland (2009)

In particular, charges $e^{A}_{i}$ and $e^{B}_{i-1}$ are identified as a new sector $e_{i}$. Moreover, individual fluxes are confined due to their non-trivial statistics with the condensed bosons, but flux pairs $m^{A}_{i}m^{B}_{i-1}$ survive the condensation as sectors labelled $m_{i}$. Sectors $e_{i}$ and $m_{i}$ have a mutual $-1$ braiding statistic, and adjacent fluxes $m_{i}$ and $m_{i+1}$ inherit the $i$ braiding statistic. Therefore, the fluxes obey the fusion rules

[TABLE]

Upon compactification, the model may be thought of as a 2D $\prod_{i=1}^{L}Z_{2}$ twisted gauge theory with type-II twists between adjacent fluxes.

III.1.2 Giant $K$-matrix

In 2D, abelian topological orders can be generically understood in terms of the $K$ matrix Chern-Simons formalism.Wen and Zee (1992) In this description, $N$ species of $U(1)$ gauge fields, $a_{I}$ with $I=1,\ldots,N$, are governed by the Lagrangian

[TABLE]

where $K$ is an $N\times N$ symmetric integer matrix, with even integers along the diagonal for bosonic systems. The quasiparticles are represented by integer vectors $l=(l_{1},\ldots,l_{N})\in\mathbb{Z}^{N}$, and have exchange statistics

[TABLE]

whereas their mutual braiding statistics are given by

[TABLE]

Quasiparticles of the form $Kl$ for $l\in\mathbb{Z}^{N}$ have trivial statistics with all other quasiparticles and thus correspond to local excitations. It is important to note that two $K$ matrices, $K$ and $K^{\prime}$, are physically equivalent if there is a unimodular matrix $W$ (i.e. with $\det W=1$) such that $K^{\prime}=W^{T}KW$. Such a transformation corresponds to a change of quasiparticle basis.

Here, we will employ the $K$ matrix formalism to describe the excitation content of the 3D condensed phase of the prior section. In particular, the structure of planons is captured by a ‘giant’ $N\times N$ $K$-matrix, whose dimension is extensive in the height of the system, and in which spatial locality of excitations in the $z$ direction is encoded in the indices of the vector $l$. In other words, the quasiparticle represented by $l=(\ldots,0,1,0,\ldots)$, with nonzero value at index $I$, is a planon constrained to move near the $xy$ plane with $z$ coordinate equal to $I$ units of the lattice spacing.

We consider the $K$ matrix with the following form in the bulk (where we have labeled the columns in the anyon basis)

[TABLE]

The inverse matrix $K^{-1}$ has the following form:

[TABLE]

The quasiparticle statistics can be read off from $K^{-1}$. Denoting by $l_{I}$ the unit vector with all entries equal to 0 except the entry at index $I$, the giant $K$ matrix corresponds precisely to the excitation content of the boson-condensed phase under the assignment $l_{2i-1}=e_{i}$ and $l_{2i}=m_{i}$. In Appendix C.3, we describe a lattice model realization of the above $K$-matrix.

III.2 Foliation structure

The foliation structure of the model can be easily understood in the $K$-matrix formalism. A single layer of $Z_{2}\times Z_{2}$ twisted gauge theory may be disentangled from the bulk via a local unitary transformation represented by the following $W$ matrix, which maps the $e_{i}$ and $m_{i}$ anyon basis to a new $\tilde{e}_{i}$ and $\tilde{m}_{i}$ basis:

[TABLE]

This $W$ matrix transforms the $K$-matrix as follows:

[TABLE]

Evidently, the transformed $K$-matrix is block diagonal. The $4\times 4$ block (for anyons $\tilde{e}^{A}$ through $\tilde{m}^{B}$), which we will call $K_{2D}$, represents a disentangled copy of 2D $Z_{2}\times Z_{2}$ twisted gauge theory. To see that this is the case, note that $K_{2D}$ has inverse

[TABLE]

On the other hand, it can easily be seen that the remaining rows and columns represent a smaller version of the original 3D model.

Because the ground state degeneracy of the system only grows with linear system size in the $z$ direction but not in the $x$ and $y$ direction, the model is 1-foliated. That is, growing the model in the $z$ direction requires the addition of 2D topological order resource layers, whereas growing the model in the $x$ or $y$ directions simply requires product state resources.

III.3 Nontrivial foliated fracton order

By examining the structure of the planon fusion group, we will demonstrate in this section that the 1-foliated model is not local unitarily equivalent to any decoupled stack of 2D topological orders, nor can it be made equivalent by adding any number of 2D topological order resource layers. In other words, the model represents a non-trivial foliated fracton phase. It is twisted in the sense that ungauging the model yields a nontrivial SSPT phase with 1 set of planar subsystem symmetries.

The situation is very similar to that of the 3-foliated model after dimension reduction. We can choose a locally complete generating set for the planons as $\{e_{i},m_{i},i=1,...,L\}$. This generating set is redundant with local redundancy relations

[TABLE]

We can start to remove the redundancy relations by eliminating the $e$’s from the generating set. However, the redundancy relations necessarily gets longer into the form $\left(m_{2}\times m_{4}\times...\times m_{2n}\right)^{2}=e_{1}\times e_{2n+1}$. Therefore, the redundancy relations cannot be locally removed and we conclude that the 1-foliated model is not equivalent to a stack of 2D layers and is hence ‘twisted’.

IV Mapping to subsystem SPT phases

The 3-foliated and 1-foliated model introduced in the previous two sections can be ‘ungauged’ into subsystem symmetry protected topological (SSPT) models You et al. (2018b); Devakul et al. (2018); Stephen et al. (2018). As the fracton models have twisted foliated fracton order, correspondingly the ungauged model has nontrivial SSPT order. In this section, we first demonstrate how the mapping works, then explain in detail our definition of SSPT order, especially a subtle difference from that given in Ref. You et al., 2018b; Devakul et al., 2018.

IV.1 The mapping

As the 3-foliated model has a ‘cage-net’ type constructionPrem et al. (2018) as discussed in section II.1, it can be ‘ungauged’ through a duality transformation similar to that described in Ref. You et al., 2018a (see also Refs. Vijay et al., 2016; Williamson, 2016). In particular, the ‘matter’ degrees of freedom $\sigma^{A}$ and $\sigma^{B}$ live at the center of the cubes in the cubic lattice. The $\sigma$s can be chosen as spin $1/2$ degrees of freedom with on-site symmetry generated by $\sigma^{A}_{x}$ and $\sigma^{B}_{x}$. Upon ‘ungauging’, the fracton Hamiltonian in Eq. 5 gets mapped to a model of the $\sigma$s with planar subsystem symmetry. The Hamiltonian is

[TABLE]

where $\tilde{O}^{A}_{c}$ and $\tilde{O}^{B}_{c}$ are obtained from $O^{A}_{c}$ and $O^{B}_{c}$ of Eq. 5 in the following way: (1) Replace the tensor product of $12$ $X^{A}$ ($X^{B}$) on the edges around the cube $c$ in $O^{A}_{c}$ ($O^{B}_{c}$) with the matter DOF $\sigma^{A}_{x,c}$ ($\sigma^{B}_{x,c}$) at the center of the cube. (2) Replace $Z^{A}_{e}$ ($Z^{B}_{e}$) on each edge with the tensor product of $4$ $\sigma^{A}_{z,c}$’s ($\sigma^{B}_{z,c}$’s) in the cubes containing the edge. Note that the phase factors in the $O^{A}$ and $O^{B}$ terms can always be expanded in the basis of $Z^{A}$ and $Z^{B}$ operators.444For example, $S=\frac{1+i}{2}+\frac{1-i}{2}Z$ and $CZ=\frac{1}{2}(1\otimes 1+Z\otimes 1+1\otimes Z-Z\otimes Z)$. Therefore, these replacement steps completely determine the $\tilde{O}$ terms from the $O$ terms. Moreover, as the $\sigma^{A}_{z}$ and $\sigma^{B}_{z}$ terms always appear as the tensor product of four around each edge, the new Hamiltonian terms are invariant under subsystem planar symmetries

[TABLE]

where $P_{xy}$, $P_{yz}$, $P_{zx}$ denote planes in the $xy$, $yz$, $zx$ direction respectively.

For the 1-foliated model, which is obtained by condensing $e^{B}_{i-1}e^{A}_{i}$ charge pairs in a stack of $Z_{2}\times Z_{2}$ twisted gauge theory models, the corresponding SSPT can be obtained from a stack of $Z_{2}\times Z_{2}$ twisted SPT Chen et al. (2013); Wang and Levin (2015) by condensing $e^{B}_{i-1}e^{A}_{i}$ charge pairs. In the SSPT model, condensing charge pairs simply means that the $Z^{B}_{2}$ symmetry of the $(i-1)$th layer is combined with $Z^{A}_{2}$ symmetry of the $i$th layer into a single symmetry generator. That is, the Hamiltonian of the system is the same as that of a decoupled stack of $Z_{2}\times Z_{2}$ twisted SPT, while the planar symmetry generators are tensor products of planar symmetry generators of the $B$ part in layer $i-1$ and the $A$ part in layer $i$.

IV.2 Definition of SSPT order

As the SSPT models are obtained by ‘ungauging’ twisted fracton models, we expect the SSPT to be ‘twisted’ as well. To be more precise, a 3D system is said to have planar subsystem symmetry protected topological (SSPT) order if

Definition 1

The model has a unique symmetric gapped ground state on any closed 3D manifold, which in the absence of symmetry can be mapped to a product state using a finite depth quantum circuit.

Two SSPT models with the same subsystem symmetry are said to have the same SSPT order if

Definition 2

The two models can be mapped to each other by adding 2D SPT layers with independent planar symmetries to each model and applying a symmetric finite depth quantum circuit.

Note that there is some subtlety in comparing the subsystem symmetry group of two models as the total symmetry group depends on system size. We consider two subsystem symmetry groups to be the same if they can be made the same by adding independent planar symmetry generators to either side.

Accordingly,

Definition 3

An SSPT model has nontrivial or ‘twisted’ SSPT order if it does not have the same SSPT order as a trivial paramagnet (a product state) with the same subsystem symmetry.

It is easy to see that once the planar symmetries are gauged, this definition of SSPT order matches the definition of foliated fracton order illustrated in Fig. 1. This definition can be generalized to models and subsystem symmetries in other dimensions in a straight forward way.

Our definition is similar but also different from that in Ref. You et al., 2018b; Devakul et al., 2018. The definition of Ref. You et al., 2018b; Devakul et al., 2018 makes use of a ‘linearly symmetric local unitary circuit’ while we use only symmetric finite depth circuits but allow the addition of SPT layers. That is, we require each unitary gate in the circuit to be symmetric while the definition in Ref. You et al., 2018b; Devakul et al., 2018 allows the individual gates to break symmetry and requires only a subsystem (linear or planar) composite of them to be symmetric. A common consequence of these two definitions is that a pure stack of lower dimensional SPTs, where the subsystem symmetry acts as a global symmetry on each of them, is considered to be a trivial SSPT. On the other hand, the ‘linearly symmetric local unitary’ equivalence is stronger. In particular, in our definition we require the added SPT to come with their own independent symmetry generators. After they are added to the total system, the total subsystem symmetry group is always enlarged. The effect of the ‘linearly symmetric local unitary’ can also be interpreted as allowing the addition of subsystem SPTs. But once added, the symmetry generator of the SPT can be identified with one of the original symmetry generators of the system, hence directly changing the SPT signature associated with that generator. Our definition of equivalence is weaker (e.g. our definition classifies more models as nontrivial) and we have chosen it so that it matches with our definition of foliated fracton order once the subsystem symmetries are gauged.

Upon gauging, the equivalence condition in Ref. You et al., 2018b; Devakul et al., 2018 is different from the foliated fracton equivalence we used in this paper. Compared to the foliated fracton equivalence, it amounts to allowing charge condensation in fracton models, because prior to gauging the symmetry group does not necessarily become larger when SPT layers are added. The 1-foliated model discussed above is a trivial SSPT phase under their definition, while it is nontrial under our definition. The 3-foliated model is likely a nontrivial SSPT phase under both definitions.

V Summary

To summarize, in this paper we demonstrate the existence of twisted foliated fracton order, i.e. 3D gapped fracton models with a foliation structure but which are inequivalent to (copies of) the X-cube model. In particular, we discussed a 3-foliated model in section II and a 1-foliated model in section III. We demonstrated the nontriviality of the models by studying the fractional excitations – the lineons and the planons – of the models. In particular, we used a dimensional reduction procedure to reduce the 3D model to a 2D model while keeping track of the locality of the planons along the reduced dimension. By studying the group structure of the local planons, we can discern the differences between stacks of 2D layers, the X-cube model and the twisted models. By using an ungauging procedure, we further mapped the twisted fracton models to nontrivial subsystem symmetry protected topological models.

Acknowledgments

We are grateful for helpful discussions with Fiona Burnell, Trithep Devakul, Yizhi You, Hao Song, Abhinav Prem, Xiuqi Ma, Tina Zhang, Roger Mong, Dominic Williamson, and Meng Cheng. W.S. and X.C. are supported by the National Science Foundation under award number DMR-1654340 and the Institute for Quantum Information and Matter at Caltech. X.C. is also supported by the Alfred P. Sloan research fellowship and the Walter Burke Institute for Theoretical Physics at Caltech. K.S. is supported by the Walter Burke Institute for Theoretical Physics at Caltech.

Appendix A The X-cube model

The X-cube model, as first discussed in Ref. Vijay et al., 2016, is defined on a cubic lattice with qubit degrees of freedom on the edges. The Hamiltonian

[TABLE]

contains two types of terms: cube terms $B_{c}$ which are products of the twelve Pauli $X$ operators around a cube $c$, and cross terms $A^{\mu}_{v}$ which are products of the four Pauli $Z$ operators at a vertex $v$ in the plane normal to the $\mu$-direction where $\mu=x,y,\text{ or }z$ (Fig. 7).

Consider an $L_{x}\times L_{y}\times L_{z}$ cubic lattice with periodic boundary conditions. The ground state degeneracy (GSD) scales linearly with the size of the system in all three directions:

[TABLE]

There are hence a large number of ‘logical operators’ that commute with all of the terms in the Hamiltonian and map one ground state to another.Slagle and Kim (2017b); He et al. (2018) An over-complete set of $X$-type logical operators is given by the set of closed string-like operators $W^{\mu}_{ij}$, which is a product of $X$ operators over all $\mu$-oriented edges with coordinates $(i,j)$ in the plane normal to $\mu$ (see Fig. 8). This set is over-complete in the sense that products of the form $W^{\mu}_{ij}W^{\mu}_{il}W^{\mu}_{kl}W^{\mu}_{kj}$ are equal to a product of some $B_{c}$ cube operators, and thus act trivially on the ground state manifold (here the four sets of coordinates lie ahbt the corners of a rectangle in the plane normal to $\mu$, as shown in Fig. 8). There are $L_{x}L_{y}+L_{y}L_{z}+L_{z}L_{x}-2L_{x}-2L_{y}-2L_{z}+3$ such relations corresponding to unique products of cube operators, thus implying Eq. (28).

Logical operators correspond to processes where particle anti-particle pairs are created out of the vacuum, wound around the torus, and then annihilated. Straight open string operators $W^{\mu}_{ij}\left(\mu_{1},\mu_{2}\right)$ anti-commute with the vertex Hamiltonian terms at the endpoints $\mu_{1}$ and $\mu_{2}$, corresponding to excitations which live on the vertices of the lattice. Here $W^{\mu}_{ij}\left(\mu_{1},\mu_{2}\right)$ is defined to be the product of $X$ operators over $\mu$-oriented edges between $\mu=\mu_{1}$ and $\mu=\mu_{2}$ with coordinate $(i,j)$ in the plane normal to $\mu$ (see Fig. 9). Conversely, acting with bent string operators introduces additional energetic costs at the corners. Therefore the particles living at the endpoints of straight open strings are energetically confined to live on a line; in this sense, they are dimension-1 particles. Vijay et al. (2016) These particles obey an unconventional fusion rule: triples of particles living along $x$-, $y$-, and $z$-oriented lines may annihilate into the vacuum. On the other hand, acting with a closed string operator around a rectangle creates an excitation at each corner of the rectangle. A pair of particles at adjacent corners may be viewed as a single dipole-like object which is itself a dimension-2 particle and is mobile in the plane normal to the edges connecting the two corners.

In addition to these string-like operators, there are membrane-like operators which are products of $Z$ operators over qubits corresponding to a membrane geometry on the dual lattice (see Fig. 9). A rectangular membrane operator anti-commutes with the cube Hamiltonian terms at its corners. A pair of adjacent corner excitations created by a rectangular membrane operator is likewise a dimension-2 dipolar particle, free to move in a plane perpendicular to its moment. A process whereby a pair of such membrane dipoles is created, separated, wound around the torus and annihilated, corresponds to a string-like $Z$-type logical operator.

Appendix B Ground state degeneracy

In this appendix, we review algorithms to compute the ground state degeneracy and entanglement entropy of a $Z_{D}$ qudit stabilizer code. Gheorghiu (2014); Fattal et al. (2004); Gottesman (1997)

Consider a stabilizer code of the form

[TABLE]

Each $s_{\alpha}$ is a product of $Z_{D}$ clock and shift operators $Z$ and $X$ where $ZX=\omega XZ$ and $\omega=e^{2\pi i/D}$. Note that $H$ is completely determined by the $k$-component integer vector $p_{\alpha}$ and $k\times 2n$ integer matrix $S$. Since we require that $H$ is a stabilizer code, any product of $s_{\alpha}$ that results in a multiple of the identity operator must be the identity operator exactly; i.e. $H$ must be frustration-free.

Multiplying one stabilizer by another or applying unitary Clifford operators to $H$ roughly corresponds to multiplying $S$ on the left or right by an invertible integer matrix, along with some additional modifications to $p_{\alpha}$. Analogous to the singular value decomposition, the Smith decomposition diagonalizes an integer matrix using invertible integer matrices. Therefore, we can compute the Smith normal form of $S$ to obtain a new integer matrix $S^{\prime}$ which is diagonal, and the Hamiltonian $H^{\prime}$ defined by $S^{\prime}$ will have the same ground state degeneracy as $H$. Since $S^{\prime}$ is diagonal, $H^{\prime}$ consists of decoupled qudits, and the ground state degeneracy of $H^{\prime}$ is trivial to calculate (and the new phases $\omega^{p^{\prime}_{\alpha}}$ do not affect the degeneracy). For the special case of $Z_{D}$ qudits with $D$ prime, the degeneracy can instead be calculated from the rank of $S$ over the field $Z_{D}$.

An algorithm to compute the entanglement of a qubit stabilizer code is discussed in Ref. Fattal et al., 2004. Similar to the ground state degeneracy calculation, the entanglement entropy is computed in terms of the rank of a matrix $S_{AB}$ over the field $Z_{D}$ when the qudit dimension $D$ is prime. For non-prime $D$, the algorithm generalizes similarly to the degeneracy calculation and the entanglement entropy is calculated from the Smith diagonals of the same matrix.

Appendix C String-membrane-net realization

The 3-foliated model in Sec. II.3.3 can also be written as a string-membrane-net (SMN) model Slagle et al. (2019). The SMN consists of two 3D $Z_{2}$ toric codes coupled to 2D $Z_{4}$ toric code (TC) layers. The coupling modifies the set of local excitations along the 2D layers, which in turn modifies the mobility of the excitations:

When a pair of charges $e_{\text{3D}}^{\text{(1)}}$ ($e_{\text{3D}}^{\text{(2)}}$) of the first (or second) 3D TC is created across a layer, a pair of charge $2e_{\text{2D}}$ (flux $2m_{\text{2D}}$) excitations is also created on the 2D TC layer. 2. 2.

When a pair of oppositely-charged $Z_{4}$ charge $\pm e_{\text{2D}}$ (or flux $\pm m_{\text{2D}}$) excitations is created on a 2D TC layer, an open $\pi$ flux string excitation of the second (first) 3D TC is also created with endpoints on the two oppositely-charged 2D excitations.

See Fig. 10 for pictures of these local excitations.

Note that the mobility of particles is determined by the set of local excitations since charges can move by creating and annihilating local excitations, such as a pair of slightly displaced excitations of opposite charge. However, exotic sets of local excitations lead to more interesting mobility rules. For example, due to the first effect above, the 3D toric code (TC) charges ($e_{\text{3D}}^{(1)}$ and $e_{\text{3D}}^{(2)}$) are fractons since they must leave behind 2D TC excitations when they pass through layers. The second effect implies that an odd number of 2D TC charges ($e_{\text{(2D)}}$) or fluxes ($m_{\text{(2D)}}$) must be attached to the endpoints of 3D TC flux strings, which implies that an odd number of 2D TC charges or fluxes are linearly confined. However, a pair of 2D TC charges (or fluxes) from two intersecting layers is a lineon because this pair is confined to the intersection of the two layers by the 3D TC flux strings.

The Hamiltonian of the string-membrane-net can be written down on very general lattices. In particular, it is possible to consider lattices where there are many qubits between the toric code layers so that one can indeed think of the Hamiltonian as 2D toric codes coupled to two 3D toric codes. In Fig. 11, we depict the simplest example where the Hamiltonian is defined on a cubic lattice in which the toric code layers are placed a single lattice spacing apart from one another.

C.1 Unitary Mapping

To show that the string-membrane-net Hamiltonian (Fig. 11) is equivalent to the cage-net Hamiltonian in Fig. 5, we will show that there is a unitary mapping between the ground spaces of the two Hamiltonians (augmented with some extra decoupled degrees of freedom).

To begin, it is convenient to replace the $Z_{2}$ qubits of the two 3D toric codes with $Z_{4}$ qudits. This will be achieved by making the following operator replacement in the string-membrane-net Hamiltonian (Fig. 11):

[TABLE]

and adding the following terms to the Hamiltonian:

[TABLE]

We have replaced the Pauli operators $\sigma^{\mu}$ and $\tau^{\mu}$ with clock and shift operators $\tilde{\sigma}^{\mu}$ and $\tilde{\tau}^{\mu}$, which have the algebra $\tilde{\sigma}^{z}\tilde{\sigma}^{x}=i\tilde{\sigma}^{x}\tilde{\sigma}^{z}$ and $\tilde{\tau}^{z}\tilde{\tau}^{x}=i\tilde{\tau}^{x}\tilde{\tau}^{z}$. The above replacement does not change the ground state since the new terms in the Hamiltonian will enforce $\tilde{\sigma}=\pm 1$ and $\tilde{\tau}=\pm 1$, and the modified Hamiltonian does not have any $\sigma^{z}$ or $\tau^{x}$ operators, only $(\tilde{\sigma}^{z})^{2}$ and $(\tilde{\tau}^{x})^{2}$ operators. Thus, its ground state is still effectively described by qubits.

The next step is to act with the unitary shown in Fig. 12, which is composed of the $Z_{4}$ controlled-X operators:

[TABLE]

The replacement in Eq. (31) and unitary in Fig. 12 map the operators of the string-membrane-net Hamiltonian (Fig. 11) to those of the cage-net Hamiltonian (Fig. 5) as follows

[TABLE]

The $A_{e}^{\text{(SMN)}}$ and $B_{p}^{\text{(SMN)}}$ operators are mapped to $\tau_{e}^{z}$ and $\sigma_{p}^{x}$. This sets $\tau_{e}^{z}=\sigma_{p}^{x}=1$ in the ground state of the new Hamiltonian. We also had to add two new terms to the Hamiltonian in Eq. (32). These new terms are mapped to

[TABLE]

But since $\tau_{e}^{z}=\sigma_{p}^{x}=1$ in the ground state, the new terms are effectively mapped to $A_{e}^{\text{(cage)}}$ and $B_{p}^{\text{(cage)}}$. Therefore, the string-membrane-net Hamiltonian (Fig. 11) and cage-net Hamiltonian (Fig. 5) both have the same ground state (up to trivial decoupled degrees of freedom).

C.2 Field theory

It is also possible to describe this model using a foliated field theory. Foliated field theories, which were introduced in Ref. Slagle et al., 2019, are field theories that explicitly couple to a foliation structure via foliation fields $e_{\mu}^{k}$.555The X-cube field theory in Ref. Slagle and Kim, 2017b was written as a foliated field theory in Ref. Slagle et al., 2019.

The Lagrangian is

[TABLE]

where $A^{k}$, $B^{k}$, $a$, and $b^{\prime}$ are 1-form gauge fields, $b$ and $a^{\prime}$ are 2-form gauge fields, $e^{k}$ are static foliation fields that describe the geometry of the foliations, and $k=1,2,..,n_{f}$ indexes the different foliation layers. The $n_{f}=3$ foliation structure of a cubic lattice is described by $e^{k}_{\mu}=\lambda\delta^{k}_{\mu}$ where $\mu=0,1,2,3$ indexes the spacetime indices and $\lambda$ is the density of foliation layers.

The Lagrangian has the following gauge invariance

[TABLE]

where $\zeta^{k}$, $\chi^{k}$, $\mu^{k}$ $\nu^{k}$, $\alpha$, and $\beta^{\prime}$ are arbitrary scalars and $\beta$ and $\alpha^{\prime}$ are arbitrary 1-forms. The Lagrangian is also self-dual under

[TABLE]

This self-duality interchanges the two 3D toric codes and interchanges the 2D toric code charge and flux sectors.

C.3 1-foliated model

In this appendix, we write down a CSS code lattice model that can describe the twisted 1-foliated K-matrix model in Eq. (20). One option would be to consider the 1-foliated version of the string-membrane-net model in Fig. 11. This appears to work, but the second toric code does not have any affect in this 1-foliated case. Thus, we will consider the simpler case of a stack of 2D $Z_{4}$ toric codes coupled to a 3D $Z_{2}$ toric code. This model is a special case of the generalized string-membrane-net model in Appendix A of Ref. Slagle et al., 2019. The model is summarized in Fig. 13.

The anyon labels in Eq. (20) have the following correspondence with the excitations of the 1-foliated string-membrane-net:

[TABLE]

The anyon $m_{2z+2}$ is equivalent to a pair of 2D fluxes on neighboring layers where the fluxes are attached to two ends of a 3D flux string. It is straightforward to check that the above anyons have the same braiding statistics as those defined in the $K^{-1}$ matrix in Eq. (21). Therefore, the lattice model in Fig. 13 is a lattice realization of the $K$-matrix in Eq. (20).