Corner Modes and Ground-State Degeneracy in Models with Gauge-Like Subsystem Symmetries

TL;DR

This paper explores models with subsystem symmetries, revealing topological ground-state degeneracy and corner modes, using an exactly-solvable lattice model and a continuum field theory approach.

Contribution

It introduces a new understanding of subsystem symmetries as gauge-like, demonstrating topology-dependent features and corner states in a novel lattice model and continuum theory.

Findings

Ground state degeneracy depends on manifold topology

Presence of localized zero-energy corner modes

Continuum field theory captures geometric effects

Abstract

Subsystem symmetries are intermediate between global and gauge symmetries. One can treat these symmetries either like global symmetries that act on subregions of a system, or gauge symmetries that act on the regions transverse to the regions acted upon by the symmetry. We show that this latter interpretation can lead to an understanding of global, topology-dependent features in systems with subsystem symmetries. We demonstrate this with an exactly-solvable lattice model constructed from a 2D system of bosons coupled to a vector field with a 1D subsystem symmetry. The model is shown to host a robust ground state degeneracy that depends on the spatial topology of the underlying manifold, and localized zero energy modes on corners of the system. A continuum field theory description of these phenomena is derived in terms of an anisotropic, modified version of the Abelian K-matrix Chern-Simons field theory. We show that this continuum description can lead to geometric-type effects such as corner states and edge states whose character depends on the orientation of the edge.

Equations

$$H=-t\sum _{\text{r}}{b}_{\text{r}+\hat{y}}^{†}{b}_{\text{r}}-\mu b_{\text{r}}^{†}{b}_{\text{r}}+h.c.,$$

$$H=$$

$$H_{1d}=$$

$$H=$$

$$A_{\text{r},\text{r}+\hat{x}}=\frac{1}{2}({\theta }_{\text{r}}-{\theta }_{\text{r}+\hat{x}})+{\alpha }_{\text{r},\text{r}+\hat{x}},$$

$$A_{\text{r},\text{r}+\hat{x}}$$

$${\theta }_{\text{r}}$$

$$H=$$

$$H=-K\sum _{\text{r}}{\sigma }_{\text{r}+\hat{y},\text{r}+\hat{y}+\hat{x}}^z{\sigma }_{\text{r},\text{r}+\hat{x}}^z.$$

$$H_{eff}=$$

$$-$$

$$H=-K\sum _{\text{r}}{\sigma }_{\text{r}}^z{\sigma }_{\text{r}+\hat{y}}^z-J\sum _{x_o}\prod _{\text{r}\in l(x_o)}{\sigma }_{\text{r}}^x{\sigma }_{\text{r}+\hat{x}}^x.$$

$$G[l(x_o)]=\prod _{\text{r}\in l(x_o)}{\sigma }_{\text{r}}^x{\sigma }_{\text{r}+\hat{x}}^x.$$

$$G[l(y_o)]=\prod _{\text{r}\in l(y_o)}{\sigma }_{\text{r}}^z{\sigma }_{\text{r}+\hat{y}}^z,$$

$$H_{compass}=-J_z\sum _{\text{r}}{\sigma }_{\text{r}}^z{\sigma }_{\text{r}+\hat{y}}^z-J_x\sum _{\text{r}}{\sigma }_{\text{r}}^x{\sigma }_{\text{r}+\hat{x}}^x.$$

$$H=-J\sum _{x_o}{\bar{\sigma }}_{x_o}^x{\bar{\sigma }}_{x_o+\hat{x}}^x.$$

$$D_s=\prod _{l(x_o)}\frac{1}{2}\left(1+\prod _{\text{r}\in l(x_o)}{\sigma }_{\text{r}}^x{\sigma }_{\text{r}+x}^x\right).$$

$$W_{l(y_o)}=\prod _{\text{r}\in l(y_o)}{\sigma }_{\text{r}}^z$$

$$W_{l(x_o)}=\prod _{\text{r}\in l(x_o)}{\sigma }_{\text{r}}^x,$$

$$H_{MF}=-\sum _{⟨\text{r}\text{r’}⟩}\sum _{i,j}{J}_{ij,\text{r}\text{r’}}{\gamma }_{\text{r}}^i{\gamma }_{\text{r’}}^j,$$

$$H=$$

$${\mathcal{L}}_E=$$

$${\mathcal{L}}_E\to {\mathcal{L}}_E+b_{\mu }{\tilde{j}}^{\mu }-\frac{in}{2\pi }{b}_{\mu }{ϵ}^{\mu \nu \lambda }{\partial }_{\nu }{a}_{\lambda }.$$

$${\mathcal{L}}_E=$$

$$\mathcal{L}=\frac{n}{2\pi }{b}_{\mu }{ϵ}^{\mu \nu \lambda }{\partial }_{\nu }{a}_{\lambda }+M^2{a}_y^2+M^2{b}_x^2.$$

$$\mathcal{S}=\frac{n}{2\pi }\int d^{3}x(b_y{\partial }_t{a}_x-a_t{\partial }_x{b}_y-b_t{\partial }_y{a}_x).$$

$$a_x$$

$$b_y$$

$$\mathcal{S}=\frac{n}{2\pi }\int dt({\bar{b}}_y{\partial }_t{\bar{a}}_x).$$

$$e^{i{\bar{b}}_y}{e}^{i{\bar{a}}_x}=e^{i\frac{2\pi }{n}}{e}^{i{\bar{a}}_x}{e}^{i{\bar{b}}_y}.$$

$$a_x={\partial }_x\varphi (x,t)$$

Full text

I Introduction

It is widely known that discrete gauge symmetries can give rise to topological order in 2+1D Kogut (1979); Lüscher (1982); Woit (1983); Kogut et al. (1983); Wen (1990); Levin and Wen (2003). This began with work on 2+1D lattice gauge theory descriptions of quantum dimer models and resonating valence Bond states Fradkin and Shenker (1979); Fradkin and Kivelson (1990); Wen (1991a); Moessner et al. (2001); Ardonne et al. (2004); Levin and Wen (2006); Fradkin (2013). Since then, there has been intense theoretical effort studying the properties of topological ordered lattice gauge theories Jersák et al. (1983); Brown et al. (1988); Stack et al. (1994); Engels et al. (1995); Greensite (2003); Wang et al. (2011). Key features of these systems include, a robust ground state degeneracy which depends on the topology of the underlying spatial manifold/lattice Wen (1989); Misguich et al. (2002); Levin and Wen (2006), fractionalized quasiparticles with unusual statistics Wen and Zee (1991); Wen (1991b, c); Kitaev (2006); Nayak et al. (2008); Levin and Gu (2012), and long-range entangled ground states Kitaev and Preskill (2006); Chen et al. (2010); Jiang et al. (2012).

A quintessential example of emergent topological order is Kitaev’s toric code model, which realizes the deconfined phase of a $\mathbb{Z}_{2}$ lattice gauge theory Kitaev (1997); Bravyi and Kitaev (1998); Kitaev (2003); Castelnovo and Chamon (2008); Tupitsyn et al. (2010). The model consists of a square lattice with spin-$1/2$ degrees of freedom defined on the links of lattice. The $\mathbb{Z}_{2}$ gauge transformation consists of flipping all spins around a single elementary plaquette. When defined on a manifold of genus $g$, the toric code system has a ground state degeneracy of $4^{g}$, which corresponds to the number of ways $\mathbb{Z}_{2}$ gauge fluxes can be threaded through non-contractible loops in the system.

Recently, there has also been significant work in understanding the role of subsystem symmetries in topological phases of matter. For a $D$ dimensional system, subsystem symmetries (also refereed to as Gauge-Like symmetries) are sets of symmetries that act independently on $d$ dimensional subregions, with $0<d<D$. Subsystem symmetries can be viewed as intermediate between gauge symmetries ([math] dimensional subregions) and global symmetries ($D$ dimensional subregions).

In connection to topology, it has been shown that subsystem symmetries can lead to unique topological phases of matter known as subsystem symmetry protected topological (SSPT) phases You et al. (2018a). SSPT phases have edge degrees of freedom that transform projectively under the subsystem symmetry. For open boundaries, SSPT’s have a subextensive ground state degeneracy protected by the subsystem symmetries. In this way SSPT’s are a subsystem generalization of (global) symmetry protected topological phases Devakul et al. (2018).

Subsystem symmetries have also been studied in connection to fractonic phases of matter Chamon (2005); Haah (2011); Nandkishore and Hermele (2018). Fracton systems are 3+1D phases of matter, characterized by immobile excitations, or excitations which are confined to sub-dimensional regions. It has been found that gauging a subsystem symmetry can lead to a fractonic phaseBravyi et al. (2011); Yoshida (2013); Williamson (2016); Vijay et al. (2016); You et al. (2018b); Vijay et al. (2016). Since fracton systems are believed to be described by rank 2 symmetric gauge theories, this field has also gained attention due to possible connections to elasticity and gravity theories Pretko (2017); Pretko and Radzihovsky (2018).

Currently, the study of subsystem symmetries has been largely based on viewing a $d$-dimensional subsystem symmetry as a global symmetry acting on $d$-dimensional subregions. However, there is also a complimentary view of a $d$-dimensional subsystem symmetry as a gauge symmetry acting on a $D-d$-dimensional subregion. For example, consider a $2d$ plane with coordinates $(x,y)$, where a subsystem symmetry acts along $1d$ $y=y_{o}$(const.) lines. Restricted to $y=y_{o}$ lines, the subsystem symmetry is a global symmetry. However, for $x=x_{o}$(const.) lines the subsystem symmetry is a local/gauge symmetry, since it only acts at the point $(x_{o},y_{o})$.

Since subsystem symmetries behave like gauge symmetries in certain subregions, we believe that salient features of lattice gauge theories may occur in systems where the low energy physics is invariant under a subsystem symmetry. In particular we ask if subsystem symmetries can lead to interesting global phenomena in the same way that gauge symmetries do in topologically ordered phases. We answer this question in the affirmative by using a $D=2$ model of bosons with a $d=1$ $U(1)$ subsystem symmetry. Using two complimentary descriptions, we show that this model has multiple ground states on a torus, which cannot be locally distinguished. Furthermore, we show that for a rectangular system with open boundaries, there are gapless degrees of freedom that are localized to the system’s corners.

This paper is organized as follows. In Section II, we construct the subsystem symmetry invariant model by using a coupled wire construction. In Section III we construct an effective projector Hamiltonian and use it to study the system. In Section IV we construct and analyze a continuum description of the subsystem symmetry invariant model. In Section V we generalize the continuum description and discuss its features. Finally, we discuss and conclude these results in Section VI.

II Subsystem Symmetry Invariant Model

To construct our subsystem symmetry invariant model, we consider an array of complex bosonic wires on a square lattice with unit directions $\hat{x}$ and $\hat{y}$. The Hamiltonian for the wire array with wires aligned parallel to the $y$-direction is given by

[TABLE]

where $b$ is a complex valued boson, and $\mu$ is a chemical potential. For a $L_{x}\times L_{y}$ lattice, this model has $L_{x}$ $U(1)$ symmetries which correspond to rotating the phase of a given wire. Formally, this symmetry operation is given by $b_{\textbf{r}}\rightarrow b_{\textbf{r}}e^{in\Lambda_{\textbf{r}}}$, where $n\in Z$, and $\Lambda_{\textbf{r}}$ is a real function that is constant along the $\hat{y}$ direction ($\Lambda_{\textbf{r}}=\Lambda_{\textbf{r}+\hat{y}}$). The factor of $n$ included in this definition is necessary for this system to have non-trivial features.

We now want to couple these wires in such a way that the $L_{x}$ $U(1)$ subsystem symmetries are preserved. To do this, we will introduce a new set of fields $A$ defined on the links that connect sites r and $\textbf{r}+\hat{x}$. These fields transform as $A_{\textbf{r},\textbf{r}+\hat{x}}\rightarrow A_{\textbf{r},\textbf{r}+\hat{x}}+(\Lambda_{\textbf{r}}-\Lambda_{\textbf{r}+\hat{x}})$ under the $U(1)$ subsystem symmetries. Introducing these fields, the Hamiltonian becomes

[TABLE]

This model now has subsystem symmetries given by $b_{\textbf{r}}\rightarrow b_{\textbf{r}}e^{in\Lambda_{\textbf{r}}}$ and $A_{\textbf{r},\textbf{r}+\hat{x}}\rightarrow A_{\textbf{r},\textbf{r}+\hat{x}}-(\Lambda_{\textbf{r}}-\Lambda_{\textbf{r}+\hat{x}})$, where $\Lambda$ is a real function that is constant along the $\hat{y}$ direction. The $t^{\prime}$ coupling in Eq. 2 can be viewed as a subsystem generalization of a gauge connection, i.e., a way of coupling the bosons such that the subsystem symmetry is preserved. This coupling has also introduced vortex configurations where the value of $A$ jumps by $2\pi/n$. The term proportional to $K$ adds an energy cost to creating these vortices. Since the $K$ terms only couple fields that are neighbors in the $\hat{y}$-direction, these vortex excitations can only propagate along the $\hat{y}$-direction.

To gain more insight into this Hamiltonian, let us restrict our attention to a line along the $\hat{x}$ direction defined as $l(y_{o})=\{\textbf{r}=(x,y)|y=y_{o}\}$ where $y_{o}$ is a constant. Let us extract the section of the Hamiltonian that acts only on $l(y_{o})$. The resulting $1d$ Hamiltonian for this subregion is

[TABLE]

where $x\equiv(y_{o},x)$. This is exactly the Hamiltonian for $1d$ charge $n$ bosons coupled to a gauge field $A$. The gauge transformations are given by $b_{x}\rightarrow b_{x}e^{in\Lambda^{\prime}_{x}}$ and $A_{x,x+1}\rightarrow A_{x,x+1}+(\Lambda^{\prime}_{x}-\Lambda^{\prime}_{x+1})$. This is exactly the subsystem transformation of the full system restricted to the $l(y_{o})$ line. So, along the $l(y_{o})$ subregion, the subsystem symmetry corresponds to a $1d$ gauge symmetry.

Motivated by this, we can consider the expectation value of the Wilson loops of the dimensionally reduced $1d$ system $W_{1d}=\exp(i\sum_{x}A_{x,x+1})$. For periodic boundary conditions, the expectation value of $W_{1d}$ can be changed by a factor of $e^{i2\pi/n}$ by threading a unit of flux through the $1d$ system. In terms of the $A$ fields, the flux threading sends $A_{x,x+1}\rightarrow A_{x,x+1}+2\pi/(nL_{x})$, for each $x$. In the full $2d$ system, $W_{1d}$ becomes the operator $W_{l(y_{o})}=\exp(i\sum_{\textbf{r}\in l(y_{o})}A_{\textbf{r},\textbf{r}+\hat{x}})$. This operator is invariant under the $U(1)$ subsystem symmetries of Eq. 1. For periodic boundaries in the $\hat{x}$ direction, we can also define a “flux insertion” operation that sends $A_{\textbf{r},\textbf{r}+\hat{x}}\rightarrow A_{\textbf{r},\textbf{r}+\hat{x}}+2\pi/(nL_{x})$ for each r. This will change the expectation value of $W_{l(y_{o})}$ by a factor of $e^{i2\pi/n}$.

It is clear that $W_{l(y_{o})}$ is similar to the Wilson loops of a $2d$ lattice gauge theory. To illustrate the similarities and differences between lattice gauge theories and Eq. 2, let us consider these systems on a torus. For a $2d$ lattice gauge theory there are two distinct non-contactable Wilson loops: one oriented in the $\hat{x}$ direction, and one oriented in the $\hat{y}$ direction. The expectation value of these loops can be changed by threading flux through the $\hat{y}$ or $\hat{x}$ directions respectively. However, for Eq. 2, the Wilson loop-like operator $W_{l(y_{o})}$ is fixed to be oriented in the $\hat{x}$ direction. As a result, the system only responds to threading flux through the $\hat{y}$ direction. Motivated by this, it will prove useful to think of Eq. 2 as a gauge theory where the Wilson loops are restricted to be oriented in the $\hat{y}$ direction, or equivalently where flux can only be inserted in the $\hat{x}$ direction.

Now let us tune $\mu$ such that there is a large boson occupancy per site. $b$ can then be replaced with the rotor variable $e^{i\theta}$, where $\theta$ corresponds to the phase of the complex boson $b$Phillips (2012). The Hamiltonian then becomes

[TABLE]

The subsystem symmetry is now given by $A_{{\textbf{r},\textbf{r}+\hat{x}}}\rightarrow A_{{\textbf{r},\textbf{r}+\hat{x}}}+(\Lambda_{\textbf{r}}-\Lambda_{\textbf{r}+\hat{x}})$, and $\theta_{\textbf{r}}\rightarrow\theta_{\textbf{r}}+n\Lambda_{\textbf{r}}$ where $\Lambda_{\textbf{r}}$ is constant along the $\hat{y}$ direction. This model is the main result of this section.

It is worth noting that due to the generalized Elitzur’s theorem Batista and Nussinov (2005), the continuous $1d$ subsystem symmetry of Eq. 4 cannot be spontaneously broken. So the ground state of Eq. 4 must be invariant under under all subsystem symmetry transformations, as must all local observables. This is similar to gauge theories, where the ground state and local observables must also be invariant under all local gauge transformations.

III Effective Projector Hamiltonian

To better study Eq. 4, it will be useful to construct an effective description in terms of an exactly solvable model of commuting projectors. The resulting model will be non-local, however it will be useful to determine key features of Eq. 4 such as ground state degeneracy, and edge physics. In Section IV, we will rederive these results using a local continuum description of Eq. 4.

We will consider the case where $t,t^{\prime}\rightarrow\infty$ while $K$ remains finite. The low energy excitations will thereby be violations of the term proportional to $K$ (vortices of $A$) in Eq. 4. To be explicit, let us consider an effective description for $n=2$. The vortices of $A$ will therefore be $\pi$-vortices, where $\exp(iA)\rightarrow-\exp(iA)$. In the large $t^{\prime}$ limit we can rewrite $A$ as

[TABLE]

where $\alpha_{\textbf{r},\textbf{r}+\hat{x}}$ is a $\pi$-valued variable ($\alpha$ only takes on values of [math] or $\pi$) that corresponds to the vortices of the $A$ field. Let us now examine how these fields transform under a subsystem symmetry transformation given by $\Lambda$ satisfying $\Lambda_{\textbf{r}}=\Lambda_{\textbf{r}+\hat{y}}$. It will be useful to decompose $\Lambda\equiv\Lambda^{s}+\Lambda^{\pi}$, where $\Lambda^{s}$ takes on values in $[0,\pi)$ and $\Lambda^{\pi}$ is a $\pi$-valued function. Under such a transformation

[TABLE]

where we have used the fact that $\theta$ is $2\pi$ periodic. Comparing Eq. 5 and 6, we see that the transformation law for $\alpha$ is $\alpha_{\textbf{r},\textbf{r}+\hat{x}}\rightarrow\alpha_{\textbf{r},\textbf{r}+\hat{x}}+(\Lambda^{\pi}_{\textbf{r}}-\Lambda^{\pi}_{\textbf{r}+\hat{x}})$. So $\alpha$ is only acted on by transformation generated by $\Lambda^{\pi}$. Since $2\Lambda^{\pi}=0\mod(2\pi)$, the transformations generated by $\Lambda^{\pi}$ form a $\mathbb{Z}_{2}$ subgroup of the full $U(1)$ group of subsystem symmetry transformations.

Because $\alpha$ is $\pi$-valued, we can identify $\exp(i\alpha)=\sigma^{z}$, where $\sigma^{z}$ is a Pauli matrix. Using Eq. 5, the Hamiltonian Eq. 4 becomes

[TABLE]

The aforementioned $\mathbb{Z}_{2}$ subsystem symmetry generated by $\Lambda^{\pi}$ flips the spins $\sigma^{z}_{\textbf{r},\textbf{r}+\hat{x}}\rightarrow-\sigma^{z}_{\textbf{r},\textbf{r}+\hat{x}}$ on an even number of columns. In terms of the spin variables, this symmetry transformation is generated by $G[l(x_{o})]=\prod_{\textbf{r}\in l(x_{o})}\sigma^{x}_{\textbf{r},\textbf{r}+\hat{x}}\sigma^{x}_{{\textbf{r}+\hat{x},\textbf{r}+2\hat{x}}}$, where $l(x_{o})=\{\textbf{r}=(x,y)|x=x_{o}\}$ (see Fig. 1).

The full Hilbert space of Eq. 7 is spanned by $\bigotimes_{\textbf{r}}\ket{\bar{\sigma}^{z}_{\textbf{r},\textbf{r}+\hat{x}}}\ket{\bar{\theta}_{\textbf{r}}}$. These are eigenstates with eigenvalues $\sigma^{z}_{\textbf{r},\textbf{r}+\hat{x}}\ket{\bar{\sigma}^{z}_{\textbf{r},\textbf{r}+\hat{x}}}=\bar{\sigma}^{z}_{\textbf{r},\textbf{r}+\hat{x}}\ket{\bar{\sigma}^{z}_{\textbf{r},\textbf{r}+\hat{x}}}$ ($\bar{\sigma}^{z}\in\pm 1$) and $\theta_{\textbf{r}}\ket{\bar{\theta}_{\textbf{r}}}=\bar{\theta}_{\textbf{r}}\ket{\bar{\theta}_{\textbf{r}}}$ ($\bar{\theta}\in[0,2\pi)$). In the $t\rightarrow\infty$ limit, we will only consider states that satisfy $\bar{\theta}_{\textbf{r}}=\bar{\theta}_{\textbf{r}+\hat{y}}$. Using this the Hamiltonian becomes

[TABLE]

In this limit the phase fluctuations are frozen out energetically and the effective model acts on the restricted Hilbert space spanned only by the spin operators $\sigma^{z}$. Formally this is a mapping that takes a state $\bigotimes_{\textbf{r}}\ket{\bar{\sigma}^{z}_{\textbf{r},\textbf{r}+\hat{x}}}\ket{\bar{\theta}_{\textbf{r}}}\rightarrow\bigotimes_{\textbf{r}}\ket{\bar{\sigma}^{z}_{\textbf{r},\textbf{r}+\hat{x}}}.$

Additionally, due to the generalized Elitzur’s theorem, all observables must be invariant under the $U(1)$ subsystem symmetries. Because of this, we should focus on just the “physical subspace” of this reduced Hilbert space, which consists of states that are invariant under the $U(1)$ subsystem symmetries generated by $\Lambda$. Under the aforementioned mapping, the physical subspace of the full Hilbert space maps to a subspace of the restricted Hilbert space that is invariant under the $\mathbb{Z}_{2}$ subsystem symmetry subgroup that acts on $\exp(i\alpha)=\sigma^{z}$. To project the restricted Hilbert space onto the corresponding physical subspace, we note that a subsystem symmetry invariant state $\ket{\psi}$ will satisfy $G[l(x_{o})]\ket{\psi}=\ket{\psi}$ for all columns $l(x_{o})$. This condition can be enforced in the low-energy subspace by adding the term $-JG[l(x_{o})]$ (with $J>0$) to the Hamiltonian Eq. 8. The resulting effective projector Hamiltonian is

[TABLE]

The low energy sector will now be invariant under the subsystem symmetry. The second term in this Hamiltonian is notably non-local. This is an artifact of projecting to the physical Hilbert space. Nevertheless, this effective model provides a simple and useful description that we can use to study the low energy features of the full system Eq. 4.

It will now be useful to simplify the lattice on which we have defined this effective spin model. Let us define a new lattice such that the sites of the new lattice are the links connecting the sites r and $\textbf{r}+{\hat{x}}$ of the original lattice. This means that the $A$ fields now live on sites instead of links. The new lattice is shown Fig. 2. After switching to the new lattice the Hamiltonian simplifies to

[TABLE]

where r are the sites on the new lattice, and $\hat{x}$ and $\hat{y}$ are now the unit directions of the new lattice. $l(x_{o})=\{\textbf{r}=(x,y)|x=x_{o}\}$ is now the set of spins along a given straight line in the $\hat{y}$ direction.

This spin model is the main result of this section. All terms in the Hamiltonian commute, and so the spin model is exactly solvable. The subsystem symmetry here is generated by

[TABLE]

This operation is shown in Fig. 3. As we can see, the non-local second term in Eq. 11 guarantees that the ground state of the system is invariant under this transformation. Eq. 11 also has a second subsystem symmetry generated by

[TABLE]

where $l(y_{o})=\{\textbf{r}=(x,y)|y=y_{o}\}$ is a line of spins in the $\hat{x}$ direction. Due to the first term in Eq. 11, the ground state will be invariant under this second subsystem symmetry as well.

Eq. 11 is similar to the quantum compass modelKugel and Khomskii (1973), a precursor to the Kitaev honeycomb modelKitaev (2003), which is given by the Hamiltonian

[TABLE]

Indeed, the quantum compass model and the spin model Eq. 11 share the same subsystem symmetries, and Eq. 11 can also arise as the effective description of the $J_{z}>J_{x}$ phase of Eq. 14 in finite sized systems. In this case, the effective $K$ will be proportional to $(J_{x}/J_{z})^{L_{y}}$. However, despite the apparent similarities, these models have different ground state properties in the thermodynamic limit. It is known that the quantum compass model has 2 phases corresponding to $J_{x}>J_{z}$ and $J_{z}>J_{x}$Dorier et al. (2005). In both phases, the number of ground states scales as $2^{L}$ for an $L\times L$ system. The $J_{x}=J_{z}$ point marks a first order phase transition that connects these two phasesChen et al. (2007). In contrast, the spin model Eq. 11 has a gapped phase with a finite number of ground states, even in the thermodynamic limit. This will be shown in the following sections.

III.1 Ground States and Excitations

The ground state of the effective spin model Eq. 11 can be found by minimizing each of the commuting terms. We can intuitively understand the nature of the ground state in the following way. The terms proportional to $K$ in Eq. 11 describe an array of decoupled Ising chains. Thus, for $J=0$, the spin model is simply an array of Ising chains in the ferromagnetic phase. In the low-energy subspace, each chain can then be characterized by a single magnetization variable $\bar{\sigma}^{z}_{x_{o}}=\langle\sigma^{z}_{\textbf{r}}\rangle_{\textbf{r}\in l(x_{o})}$.

The terms proportional to $J$ in Eq. Eq. 11 flip all spins on a pair of the neighboring Ising chains (see Fig. 3), i.e., each term flips a pair of magnetizations, e.g., $\bar{\sigma}^{z}_{x_{o}}$ and $\bar{\sigma}^{z}_{x_{o}+\hat{x}}$. Let us define the operator $\bar{\sigma}^{x}_{x_{o}}=\prod_{\textbf{r}\in l(x_{o})}\sigma^{x}$. Since $\bar{\sigma}^{z}_{x_{o}}=\langle\sigma^{z}_{\textbf{r}}\rangle_{\textbf{r}\in l(x_{o})}$, $\bar{\sigma}^{x}_{x_{o}}\bar{\sigma}^{z}_{x_{o}}=-\bar{\sigma}^{x}_{x_{o}}\bar{\sigma}^{z}_{x_{o}}$. In terms of $\bar{\sigma}$, the Hamiltonian Eq. 11 becomes

[TABLE]

This Hamiltonian is just another ferromagnetic Ising chain, with the ferromagnetism oriented in the $x$-direction. So the effect of the term proportional to $J$ in Eq. 11 is to orient the magnetization of the original Ising chains. In particular, if we start with a ground state for $J=0$, we can determine the ground state for $J>0$ by acting on the $J=0$ ground state with the operator

[TABLE]

To see this, consider a state $\ket{\psi}$ that minimizes Eq. 11 with $J=0$. Then $\sigma^{z}_{\textbf{r}}\sigma^{z}_{{\textbf{r}}+y}\ket{\psi}=\ket{\psi}$ for all r. Since $\sigma^{z}_{\textbf{r}}\sigma^{z}_{{\textbf{r}}+y}D_{s}\ket{\psi}=D_{s}\sigma^{z}_{\textbf{r}}\sigma^{z}_{{\textbf{r}}+y}\ket{\psi}=D_{s}\ket{\psi}$, $D_{s}\ket{\psi}$ minimizes all $K$ terms in Eq. 11. It is also true that $(\prod_{\textbf{r}\in l(x_{o})}\sigma^{x}_{\textbf{r}}\sigma^{x}_{{\textbf{r}}+\hat{x}})D_{s}=D_{s}$, for all $x_{o},$ and by extension, $(\prod_{\textbf{r}\in l(x_{o})}\sigma^{x}_{\textbf{r}}\sigma^{x}_{{\textbf{r}}+\hat{x}})D_{s}\ket{\psi}=D_{s}\ket{\psi}$. So $D_{s}\ket{\psi}$ also minimizes all $J$ terms in Eq. 11. $D_{s}\ket{\psi}$ thereby minimizes the entire Hamiltonian with $J>0,$ and is the ground state.

We note here that $D_{s}$ is in fact exactly the projection operator that projects the restricted Hilbert space of Eq. 8 to the subsystem symmetry invariant physical subspace of the restricted Hilbert space. As we shall demonstrate below, the number of ground states will depend on the topology of the lattice. The excited states of the spin model are characterized by having either $\sigma^{z}_{\textbf{r}}\sigma^{z}_{{\textbf{r}}+\hat{y}}=-1$ or $\prod_{\textbf{r}\in l(x_{o})}\sigma^{x}_{\textbf{r}}\sigma^{x}_{\textbf{r}+\hat{x}}=-1$, which have an excitation energy of $2K$ and $2J$ respectively.

III.2 Ground State Degeneracy

A key feature of the subsystem symmetry invariant model Eq. 4 is the existence of multiple ground states that cannot be locally distinguished. We will demonstrate this by considering the effective spin model Eq. 11 on a torus. To find the number of ground states, we will identify operators that commute with the Hamiltonian, and use them to label the degenerate ground states. The non-trivial operators that commute with the Hamiltonian Eq. 11 are

[TABLE]

where $l(y_{o})=\{\textbf{r}=(x,y)|y=y_{o}\}$ is a closed loop in the $\hat{x}$ direction, and $l(x_{o})=\{\textbf{r}=(x,y)|x=x_{o}\}$ is a closed loop in the $\hat{y}$ direction. On a torus, the $l(y_{o})$ and $l(x_{o})$ loops will be the two cycles that define the torus. These loops are shown in Fig. 4. For an $L_{x}\times L_{y}$ torus, the total number of $W_{l(y_{o})}$ operators is $L_{y}$ and the number of $W_{l(x_{o})}$ operators is $L_{x}$. Since $W^{2}_{l(x_{o})}=W^{2}_{l(y_{o})}=1$, both loop operators are $\mathbb{Z}_{2}$ operators. We can identify the symmetry operator $G[l(x_{o})]$ as the product of the neighboring loop operators $W_{l(x_{o})}$ and $W_{l(x_{o}+\hat{x})}$, and similarly identify $G[l(y_{o})]$ as the product of $W_{l(y_{o})}$ and $W_{l(y_{o}+\hat{y})}$.

All loops $l(x_{o})$ and $l(y_{o})$ on the torus intersect once, and so all $W_{l(y_{o})}$ operators anti-commute with all $W_{l(x_{o})}$ operators. The minimum dimension needed to represent this anti-commuting algebra is 2, leading to 2 distinct ground states. If we were to diagonalize the ground state subspace to label them by their $W_{l(y_{o})}$ eigenvalue, then the 2 ground states would be related by the acting on a given ground state with an operator $W_{l(x_{o})}$. Since $W_{l(y_{o})}$ and $W_{l(x_{o})}$ are both non-local operators, this degeneracy is robust to local perturbations.

The degeneracy can also be found by counting the number of constraints for $L_{x}\times L_{y}$ spins on a torus. Let us first consider the terms proportional to $K$ in Eq. 11. These terms describe a system of $L_{x}$ Ising chains with periodic boundaries. Each chain contributes $L_{y}-1$ unique constraints, leading to $L_{x}(L_{y}-1)$ unique constraints from the $K$ terms in Eq. 11. The terms proportional to $J$ in Eq. 11 then give $L_{x}-1$ unique constraints. Since all terms in Eq. 11 commute, all these constraints can be simultaneously satisfied, leading to $L_{x}(L_{y}-1)+L_{x}-1=L_{x}\times L_{y}-1$ constraints in total. There is thereby $1$ net free spin degree of freedom which corresponds to the $2$ ground states that were previously identified.

It is useful to compare these results to the case of a $\mathbb{Z}_{2}$ lattice gauge theory on a torus. In $\mathbb{Z}_{2}$ lattice gauge theory models, there are 2 additional ground states on a torus (for a total of 4 ground states)Kitaev (2006). These 2 additional ground states occur since non-contractible loops of $\sigma^{x}$ operators oriented in the $\hat{x}$ direction, and non-contractible loops of $\sigma^{z}$ oriented in the $\hat{y}$ direction also commute with the $\mathbb{Z}_{2}$ lattice gauge theory Hamiltonian, and anti-commute with each other. These operators do not commute with the spin model Eq. 11, and so the number of ground states is reduced to 2.

On a sphere all string operators $W_{l(y_{o})}$ and $W_{l(x_{o})}$ commute, and so the ground state of Eq. 11 on a sphere is unique. We also show this explicitly in Appendix A by counting constraints. This topology-dependent degeneracy is reminiscent of the topological ground state degeneracy found in topological ordered systems, though it is important to note that our spin model has a non-local constraint. The non-locality will be removed when we discuss the continuum limit.

III.3 Open Boundaries and Corner Modes

We shall now consider the system with open boundaries. For simplicity, we shall take the lattice to be an $L_{x}\times L_{y}$ rectangle with open boundaries. For this geometry, the terms proportional to $K$ in Eq. 11 give $L_{x}(L_{y}-1)$ constraints, and the terms proportional to $J$ give $L_{x}-1$ constraints, leading to $L_{x}\times L_{y}-1$ constraints which can be simultaneously satisfied. There is then a single free spin 1/2 degree of freedom, leading to 2 ground states.

In the string picture, this can be seen by the anti-commutation between the zero energy operators $W_{l(y_{o})}$ and $W_{l(x_{o})}$ from Eq. 17 where $l(y_{o})$ (resp. $l(x_{o})$) is now a string in the $\hat{x}$ (resp. $\hat{y}$) direction stretching from one boundary to the other. Since the system has open boundaries, the string operators do not have to form closed loops to commute with the Hamiltonian and be invariant under the subsystem symmetries of the model. Since all $W_{l(y_{o})}$ operators anti-commute with all $W_{l(x_{o})}$ operators there are degenerate 2 ground states.

Furthermore, these 2 ground states correspond to anti-commuting corner degrees of freedom. To show this, we will switch to a Majorana representation of the spin-1/2 degrees of freedomKitaev (2003). This is done by introducing 4 Majorana degrees of freedom at each site $\gamma^{1}$, $\gamma^{2}$, $\gamma^{3}$, and $\gamma^{4}$. The spin degrees of freedom then become $\sigma^{x}=i\gamma^{1}\gamma^{2}$, $\sigma^{y}=i\gamma^{1}\gamma^{3}$, and $\sigma^{z}=i\gamma^{1}\gamma^{4}$, with the local constraint that $\gamma^{1}\gamma^{2}\gamma^{3}\gamma^{4}=1$. Setting $K=J$, the spin model can be described by the mean field Hamiltonian

[TABLE]

where $J_{ij,\textbf{r}\textbf{r'}}$ corresponds to the Majorana dimerization pattern given in Fig 5. The construction of this mean field Hamiltonian is outlined in Appendix B. The ground state of the spin model can then be found by projecting the mean field ground state of this Majorana Hamiltonian onto the physical states using the projector $P=\sum_{\textbf{r}}\frac{1}{2}(1+\gamma^{1}_{\textbf{r}}\gamma^{2}_{\textbf{r}}\gamma^{3}_{\textbf{r}}\gamma^{4}_{\textbf{r}})$, and then using the above identification between the spin operators and Majorana fermions.

For the rectangular geometry, there are 4 free Majorana degrees of freedom located at the corners of the lattice (see Fig. 5). This leads to 4 ground states, but only 2 are physical after projecting onto the physical states. The reduction in ground states can be viewed as a consequence of the full model being bosonic, i.e., fermion parity even. For the rectangular geometry, the degrees of freedom for the spin model are thereby zero energy corner operators leading to a robust ground state degeneracy. We note that this particular dimerization pattern and corner mode configuration is similar to an insulator model presented in Ref. Benalcazar et al., 2017a that has corner charge, but vanishing quadrupole moment. Indeed, the corner modes considered here are similar to what is found in higher order topological insulatorsBenalcazar et al. (2017b); Langbehn et al. (2017); Schindler et al. (2018); Khalaf (2018). However, for the spin model Eq. 11, the corner operators are non-local. This is because individual Majorana corner operators do not commute with the projector operator $P$. Only pairs of Majorana corner operators commute with $P$ and are physical.

IV Continuum Theory

We now seek a complementary continuum description of Eq. 4. First, we note that Eq. 4 is the low energy description of

[TABLE]

where $A$ fields are now defined along both $\hat{x}$ and $\hat{y}$ oriented links, and $\langle\textbf{r}\textbf{r'}\rangle$ are neighboring sites. The sum over $p$ is over plaquettes with corners $i,j,k,l,$ and $F_{p}=A_{i,j}-{A_{k,l}}+A_{i,k}-A_{j,l}$. In the low energy ($\bar{m}^{2}\rightarrow\infty$) limit, $A_{\textbf{r},\textbf{r}+\hat{y}}$ is pinned to be [math] by the cosine term. Upon substituting this into Eq. 19, the Hamiltonian reduces to Eq. 4 with $t^{\prime}=t$. This model has the subsystem symmetry given by $A_{\textbf{r},\textbf{r'}}\rightarrow A_{\textbf{r},\textbf{r'}}+(\Lambda_{\textbf{r}}-\Lambda_{\textbf{r'}})$, and $\theta_{\textbf{r}}\rightarrow\theta_{\textbf{r}}+n\Lambda_{\textbf{r}}$ where $\Lambda_{\textbf{r}}=\Lambda_{\textbf{r}+\hat{y}}$. This is the same symmetry as in Eq. 4. We note that $A_{\textbf{r},\textbf{r}+\hat{y}}$ is invariant under these transformations. Eq. 19 is the Hamiltonian for bosons minimally coupled to a vector field $A$ with an additional mass term for the fields $A$ oriented along the $y$ direction. It is worth explicitly stating that this model is not a gauge theory due to the additional mass term for $A_{\textbf{r},\textbf{r}+\hat{y}}$.

The continuum description of Eq. 19 in Euclidean space is

[TABLE]

where $\rho$ is a constant, and we have included a current $j^{\mu}$ that couples to the fluctuations of the $A$ field. To study the dynamics of the phase $\theta$, we will introduce the variable $a_{\mu}=-\frac{1}{n}\partial_{\mu}\theta$, and shift $A_{\mu}\rightarrow A_{\mu}-a_{\mu}$. After this Eq. 20 becomes

[TABLE]

where we have introduced $M^{2}\equiv\bar{m}^{2}+\rho n^{2}$. Since the $A$ field is now massive, it can be integrated out, leaving a theory just in terms of $a$. The $a$ field is now also coupled to $j$, causing the excitations of the $a$ field to have a corresponding current $j$. As desired, this model has a subsystem symmetry where $a_{\mu}\rightarrow a_{\mu}+\partial_{\mu}\Lambda_{a}(x,t)$ and $\Lambda_{a}$ is function of $x$ and $t$ only.

From the equations of motion for $a_{y}$, we have that $a_{y}\propto\frac{1}{M^{2}}$. As a result, in the low energy limit where $M^{2}\rightarrow\infty$, $a_{y}$ vanishes. In this limit, $j^{y}$ (the current in the $y$ direction) is removed from the theory, and the only current is in the $x$ direction ($j^{x}$). This means that the excitations of the $a$ field will only move in the $x$ direction. This is a form of sub-dimensional dynamics, where the excitations are only able to move in certain lower dimensional subregions.

We will also need to consider the vortex dynamics of the $a$ field. This is done by introducing a vortex current $\tilde{j}^{\mu}$, and setting $\tilde{j}^{\mu}=\frac{in}{2\pi}\epsilon^{\mu\nu\lambda}\partial_{\nu}a_{\lambda}$. To enforce this constraint, we will introduce the field $b_{\mu}$ as a Lagrange multiplier for Eq. 20

[TABLE]

In this construction, there are vortex currents in both the $x$ and $y$ directions ($\tilde{j}^{x}$ and $\tilde{j}^{y}$). However, in the original lattice model Eq. 4, the vortices were only able to move in the $y$ direction. To remedy this, we will add the term $M^{2}b^{2}_{x}$. The equation of motion for $b_{x}$ then gives that $b_{x}\propto\frac{1}{M^{2}}$, and in the low energy limit ($M^{2}\rightarrow\infty$) $b_{x}\rightarrow 0$. In this limit, the $\tilde{j}^{x}$ vortex current is removed from the theory, and there is only a vortex current in the $y$-direction ($\tilde{j}^{y}$). As a result, the vortices of the $a$ field are confined to move in the $y$-direction as in the lattice model.

After adding the $b$ field, integrating out the massive $A$ field, and keeping only the long wavelength contributions, the Lagrangian density becomes

[TABLE]

This model is the main result of this section. It is worth noting that this theory has the form of a mutual Chern-Simons theory with additional mass terms for $a_{y}$ and $b_{x}$ Bardeen and Stephen (1965); Hansson et al. (2004); Diamantini et al. (2006). This observation will be allow us to generalize this model in Section V.

In Eq. 23, it is also apparent that there is a second subsystem symmetry where $b_{\mu}\rightarrow b_{\mu}+\partial_{\mu}\Lambda_{b}(y,t)$ and $\Lambda_{b}$ is only a function of $y$ and $t$. This is the same as the subsystem symmetry generated by $G[l(y_{o})]$ (Eq. 13) in the effective projector Hamiltonian.

IV.1 Ground State Degeneracy

We will now calculate the ground state degeneracy of the continuum model on a torus. To do this, we will first rotate back to Minkowski space, and set the currents $j=\tilde{j}=0$,

[TABLE]

From the equations of motion for $a_{y}$ and $b_{x}$, we have that $a_{y}\propto b_{x}\propto\frac{1}{M^{2}}$. At low energies, $a_{y}\rightarrow 0$ and $b_{x}\rightarrow 0$, and the action becomes

[TABLE]

If we minimize the action with respect to $b_{t}$ and $a_{t}$ we find the equations of motion $\partial_{y}a_{x}=0$ and $\partial_{x}b_{y}=0$. On a torus, these equations of motion are solved by

[TABLE]

Here $\phi$ is a function of $x$ and $t$ only and is periodic on the torus, $\theta$ is a function of $y$ and $t$ and is periodic on the torus, and $\bar{a}_{x}$ and $\bar{b}_{y}$ are functions of $t$ only. $L_{x,y}$ are the length dimensions of the torus.

After substituting these terms into Eq. 25 and integrating over the $x$ and $y$ coordinates, the action reduces to

[TABLE]

Using canonical commutation relations, we have that $[\bar{b}_{y},\bar{a}_{x}]=i2\pi/n$. Since $\bar{b}_{y}$ and $\bar{a}_{x}$ are $2\pi$ periodic variables, the observables are $\exp(i\bar{b}_{y})$ and $\exp(i\bar{a}_{x})$, which obey the commutation relationship,

[TABLE]

In order to satisfy this operator algebra, there must be $n$ ground states. This is consistent with what was found using the effective projector model with $n=2$. We note that for a conventional mutual Chern-Simons theory the ground state degeneracy would be $n^{2}.$

IV.2 Corner Modes

To find the edge degrees of freedom for a system with open boundaries we will use the low energy description with no external currents in Minkowski space (Eq. 25). For a rectangular system with open boundaries, the equations of motion for $a_{x}$ and $b_{y}$ are solved by

[TABLE]

Using this, the action becomes

[TABLE]

which is a total derivative for both $x$ and $y$. If the system is defined on the rectangle $x_{0}\leq x\leq x_{1}$ and $y_{0}\leq y\leq y_{1}$, the action becomes

[TABLE]

This action describes localized operators $\eta_{j,k}\equiv\exp(i\phi(x_{j},t)-i\theta(y_{k},t))$, which are defined at the corners of the system $(x_{j},y_{k})$. Since the Hamiltonian corresponding to the action Eq. LABEL:eq:CornerAct vanishes, the $\eta_{i,j}$ are zero energy operators. It should be noted that there is a redundancy in the corner mode description, since $\eta_{0,0}\eta_{1,1}\eta^{\dagger}_{1,0}\eta^{\dagger}_{0,1}=1$.

Using the canonical commutation relationships from Eq. LABEL:eq:CornerAct, the $\eta_{j,k}$ operators satisfy the algebra

[TABLE]

Naively this would lead to $n^{2}$ ground states. However, if the constraint $\eta_{0,0}\eta_{1,1}\eta^{\dagger}_{1,0}\eta^{\dagger}_{0,1}=1$ is taken in account there are actually only $n$ ground states. This agrees with what was found in using the effective model for $n=2$. As opposed to $2+1D$ Abelian Chern-Simons field theories, where the edge theory is a $1+1D$ CFT Balachandran et al. (1991); Dunne (1999); Wen (2004); Fujita et al. (2009); Lu and Vishwanath (2012), the edge theory of the subsystem symmetry invariant model is given by $0+1D$ corner modes.

V Generalized Continuum Theory

To generalize the continuum description to include more vector fields, we note that Eq. 23 has the form of a Chern-Simons field theory with $K$-matrix $2\sigma^{x},$ and mass terms for $a_{y}$ and $b_{x}$. Using this observation, we can generalize the continuum description of the subsystem symmetry invariant system by using the Lagrangian

[TABLE]

where $K$ is a $D\times D$ symmetric integer valued matrix. We will take $M_{x,y}$ to be diagonal with all entries either $m$ or [math]. As we shall see, in order for the canonical quantization to be consistent, we shall require the number of zero eigenvalues of $M_{x}$ and $M_{y}$ to be equal, i.e., dim(ker($M_{x}$)) = dim(ker($M_{x}$)) $=\mathcal{k}\leq D$.

Minimizing the action with respect to $a^{I}_{x}$ and $a^{I}_{y}$ gives the equations of motion

[TABLE]

where $F_{\mu\nu}^{J}=\epsilon^{\mu\nu}\partial_{\mu}a_{\nu}^{J}$. Let us write $M_{x,y}^{IJ}=m\times\bar{M}_{x,y}^{IJ}$, where $\bar{M}_{x,y}^{IJ}$ is a diagonal matrix of $1$’s and [math]’s. Using $\bar{M}$, the equations of motion are

[TABLE]

At low energies $m\rightarrow\infty$, and

[TABLE]

So all vector fields $a^{J}_{x,y}$ not in the respective kernels of $M_{x,y}$ are set to [math] as $m\rightarrow\infty$. Setting $m\rightarrow\infty$ is thereby equivalent to projecting $a^{J}_{x,y}$ on to the respective $\mathcal{k}$ dimensional kernels of $M_{x,y}$. Since $M_{x,y}$ is a diagonal matrix, the kernel is spanned by a set of $\mathcal{k}$ = dim(ker($M_{x,y}$)) unit vectors. This means we can project onto the kernels of $M_{x,y}$ with $\mathcal{k}\times D$ matrices $V_{x,y}$, the rows of which are the unit vectors that span the kernels of $M_{x,y}$.

The theory with $m\rightarrow\infty$ can then be expressed as follows. Define the reduced $K$ matrices as

[TABLE]

and the vectors

[TABLE]

The effective Lagrangian in the $m\rightarrow\infty$ limit is then

[TABLE]

where we have explicitly included the sum for clarity.

V.1 Quantization

In order to consistently, canonically quantize Eq. 41, we need the following equations to be consistent

[TABLE]

To simplify this, we will define the matrix $A^{kj}=[\tilde{a}_{x}^{k},\tilde{a}_{y}^{j}]$. Eq. 42 then becomes (using $K_{yx}^{ij}=K_{xy}^{ji}$)

[TABLE]

Let us consider the case where dim(ker($M_{x}$)) = $\mathcal{k}$ and dim(ker($M_{y}$)) = $\mathcal{k}^{\prime}$. Then $K_{xy}=K^{T}_{yx}$ is a $\mathcal{k}\times\mathcal{k}^{\prime}$ matrix, and $A^{kj}$ is a $\mathcal{k}^{\prime}\times\mathcal{k}$ matrix. Summing over the $i$ and $k$ indices in Eq. 43 gives

[TABLE]

Since $Tr(K_{xy}A)=Tr(AK_{xy})$, $\mathcal{k}=\mathcal{k}^{\prime}$ in order for the quantization conditions to be consistent. This confirms our earlier assertion that we must have dim(ker($M_{x}$)) = dim(ker($M_{y}$)).

If $\det(K_{xy})\neq 0$, Eq. 42 is solved by $[\tilde{a}_{x}^{i},\tilde{a}_{y}^{j}]=2\pi i(K_{xy}^{-1})^{ij}$. If $\det(K_{xy})=0$, the inverse of $K_{xy}$ will not be well defined, and the commutation relations will be ambiguous. Because of this, we shall assume that $\det(K_{xy})\neq 0$ from here on. It is worth noting that $K_{xy}$ must be square, but $K_{ty}$ and $K_{tx}$ do not need to be square.

V.2 Ground State Degeneracy

We will now show that the ground state degeneracy on a torus is $|\det(K_{xy})|$ (which is valid because $K_{xy}$ is a square matrix). Let us minimize the action by setting the functional derivative of Eq. 41 with respect to $\tilde{a}_{t}^{i}$ equal to [math]. The resulting equations of motion are

[TABLE]

Since we are only concerned with global features of the system, we will use the solutions $\tilde{a}^{i}_{x,y}=\bar{a}^{i}_{x,y}/L_{x,y}$ where $\bar{a}^{i}_{x,y}$ is only a function of $t$, and $L_{x,y}$ are the lengths of the torus. Other solutions represent local fluctuations and do not contribute to global features.

Plugging these solutions into Eq. 41 and integrating over the $x$ and $y$ coordinates, we arrive at the action

[TABLE]

From this we can calculate the algebra satisfied by the observables $\exp(i\bar{a}^{i}_{x,y})$. Using that $[\bar{a}^{i}_{x,y},\bar{a}^{j}_{x,y}]=2\pi i(K_{xy}^{-1})^{ij}$, the minimum dimension needed to satisfy the algebra of the $\exp(i\bar{a}^{i}_{x,y})$ operators is $|\det(K_{xy})|$. This leads to a ground state degeneracy of $|\det(K_{xy})|$ on a torus.

V.3 Edge and Corner States and Example

We will illustrate some of the ground state degeneracy and edge state possibilities using an example case. Consider a $4\times 4$ $K$-matrix, and mass matrices $M_{x},M_{y}$ having kernel dimension equal to $2$. For the first example we will choose the $K$-matrix and mass matrices

[TABLE]

The corresponding fields will be labeled as $a^{i}_{\mu}$ with $i=1,2,3,4$ and $\mu=(t,x,y)$. The $2\times 4$ $V_{x,y}$ matrices for the theory are

[TABLE]

The reduced $K$ matrices are given by

[TABLE]

In the low energy limit, the Lagrangian density is given by

[TABLE]

The canonical quantization commutation relations are $[a^{1}_{y},a^{2}_{x}]=i\pi$ and $[a^{4}_{y},a^{1}_{x}]=i2\pi$. Minimizing the action with respect to $a^{i}_{t}$, we determine the equations of motion

[TABLE]

Let us consider the case where the system is put on a torus with side lengths $L_{x}$ and $L_{y}$. If we ignore local fluctuations of the fields, the equations of motion can be solved by $a^{1}_{x}=\bar{a}^{1}_{x}(t)/L_{x}$, $a^{1}_{y}=\bar{a}^{1}_{y}(t)/L_{y}$, $a^{2}_{x}=\bar{a}^{2}_{x}(t)/L_{x}$, and $a^{4}_{y}=\bar{a}^{4}_{y}(t)/L_{y}$. Since $\bar{a}^{i}_{\mu}$ is a $2\pi$ periodic variable, we will consider the operators $e^{i\bar{a}^{i}_{\mu}}.$ Using the canonical commutation relations, $e^{i\bar{a}^{y}_{x}}$ and $e^{i\bar{a}^{2}_{x}}$ anti-commute, while all other terms commute. This means that there will be 2 ground states. This agrees with $|\det(K_{xy})|=2$.

Let us now consider the edge and corner modes of this system. The equations of motion Eq. 51 are solved by $a^{1}_{x,y}=\partial_{x,y}\phi^{1}(x,y,t)$, $a^{2}_{x}=\partial_{x}\phi^{2}(x,t)$, and $a^{4}_{y}=\partial_{y}\phi^{4}(y,t)$. Plugging these solutions into Eq. 50, the action becomes

[TABLE]

To illustrate the edge modes of this system, it will be useful to consider a half plane $x\leq 0$. If we assume that the fields vanish at spatial infinity, we can rewrite the action as

[TABLE]

This describes a non-chiral boson propagating along the $x=0$ edge. If we instead consider the $y\leq 0$ half plane the action is

[TABLE]

This describes a different non-chiral boson that propagates along the $y=0$ edge. To see that this is in fact a different non-chiral boson, we note that Eq. 53 and 54 describe a $U(1)_{1}$ and $U(1)_{2}$ non-chiral boson CFT respectively.

It will also be useful to consider a quarter plane geometry $x\leq 0$ and $y\leq 0$. In this geometry, the action becomes

[TABLE]

As expected, the first line of Eq. 55 describes a non-chiral boson propagating along the $x=0,y\leq 0$ boundary, while the second line describes a different type of non-chiral boson propagating along the $x\leq 0,y=0$ boundary. Since $\phi^{1}(0,y,t)$ and $\phi^{1}(x,0,t)$ coincide at the point $(x,y)=(0,0)$, the two non-chiral bosons are coupled at this point. This means that there will be a partially-transmitive boundary connecting the two edges.

Let us now consider a different example using the same $K$ matrix as in Eq. 47, but with the mass matrices

[TABLE]

Following the same analysis as before,

[TABLE]

This gives a ground state degeneracy of $|$det$(K_{yx})|=8$.

The equations of motion for these mass matrices are

[TABLE]

The equations of motion are solved by $a^{1}_{y}=\partial_{y}\phi^{1}(y,t)$, $a^{2}_{x}=\partial_{x}\phi^{2}(x,t)$, and $a^{3}_{x,y}=\partial_{x,y}\phi^{3}(x,y,t)$.

For a half plane with $x\leq 0$, the edge action becomes

[TABLE]

This edge describes a $U(1)_{4}$ chiral boson. For a half plane with $y\leq 0$, the edge action is similarly

[TABLE]

which describes the same $U(1)_{4}$ chiral boson. Finally, for the quarter plane $x\leq 0$, $y\leq 0$, the action is

[TABLE]

The first two parts of the edge action Eq. 59 describe the previously shown chiral edge modes. The final term describes a zero energy excitation $\exp(i[\phi^{1}(0,t)-\phi^{2}(0,t)])$ which is bound to the $(x,y)=(0,0)$ corner of the system. The boundary of this system thereby has coexisting propagating chiral edge modes, and localized corner modes.

For a general theory given by Eq. 33, there will be corner modes if there exists $i$ and $j$ ($i\neq j$) such that $K^{ij}\neq 0$, $M^{ii}_{x}=m$, $M^{jj}_{y}=m$ $M^{jj}_{x}=0$, and $M^{ii}_{y}=0$. Setting $a^{i}_{y}=\partial_{y}\phi^{i}(y,t)$, and $a^{j}_{x}=\partial_{x}\phi^{j}(x,t)$, the localized mode for a corner at ($x_{o},y_{o}$) is given by $\exp(i[\phi^{i}(y_{o},t)-\phi^{j}(x_{o},t)])$.

These examples highlight the unusual edge/boundary physics that can occur in models of the form Eq. 33. In general, the edge theory for a given $K$ matrix and pair of mass matrices will consists of both propagating edge modes, and localized corner modes. Furthermore, for a given system, the edge theory may depend on the orientation of the edge since the theory can naturally support anisotropy.

VI Discussion and Conclusion: Gauging Subsystem Symmetry and Open Problems

In this work we have shown that invariance under a subsystem symmetry can lead to a topology-dependent ground state degeneracy, and corner modes. We established this using both a exactly-solvable, but non-local, spin model, and a continuum field theory description. From this, we have shown that global, topology-dependent features can exist beyond the established paradigm of gauge symmetries.

In recent literature, it has been shown that gauging subsystem symmetries of certain models can lead to fractonic phases of matterBravyi et al. (2011); Yoshida (2013); Williamson (2016); Vijay et al. (2016); You et al. (2018b); Vijay et al. (2016). For the model considered here, gauging the subsystem symmetry in our non-local spin model leads to a (local) $\mathbb{Z}_{2}$ quantum double model, i.e., the toric code. In Appendix C, we explicitly gauge the subsystem symmetry of the effective projector Hamiltonian Eq. 11, and show that this exactly leads to the toric code Hamiltonian. We can also see this in the continuum by setting $M^{2}=0$ in Eq. 23. Because of this, the subsystem symmetry invariant model we constructed can be thought of as a $\mathbb{Z}_{2}$ quantum double model, where we have “un-gauged” the $\mathbb{Z}_{2}$ gauge symmetry along a certain direction, and reduced it to a $\mathbb{Z}_{2}$ subsystem symmetry.

It still remains to be seen how this construction generalizes to higher dimensions and different subsystem symmetries. In particular if there are general topological features of a $D$ dimensional system with a $d<D$ dimensional subsystem symmetry. The Mermin-Wagner theorem would seem to constrain $d<2$ for local theories, but exact details still need to be determined.

Additionally, it is unknown what kind of classifications exist for these systems. It is known that topologically ordered systems can be classified based on their modular $S$ and $T$ matricesLevin and Wen (2005); Chen et al. (2010); Gu et al. (2015). Since the systems discussed here are not invariant under modular transformations, we cannot define these $S$ and $T$ matrices in an analogous way. Currently, to our knowledge, there is no structure which performs to same role for the subsystem invariant systems, and so classification of these systems remains an open question.

Acknowledgements

We would like to thank H. Goldman and M. Lin for useful conversations. TLH acknowledges support from the US National Science Foundation under grant DMR 1351895-CAR.

Appendix A Ground State of Effective Projector Model of a Sphere

On a manifold with the topology of a sphere (genus 0), the effective projector Hamiltonian Eq. 11 has a unique ground state. We can show this explicitly by giving the lattice the topology of a sphere. To do this we will take two commuting copies of the system stacked on top of each other. The two copies are made into a sphere by ”sewing” the copies together at the edges as shown in Fig. 6. After this, the Hamiltonian becomes

[TABLE]

where the $\pm$ indexes the two stacked copies of the system, $y{-(+)}$ are the sites on the bottom (top) edge and $x{-(+)}$ are sites on the right (left) edge. As before, all terms present in the Hamiltonian commute. The ground state is thereby determined by minimizing each term individually.

Let us now count the constraints for this system. Consider a sphere created from sewing together two $L\times L$ lattices, leading to $2L^{2}$ spins. The first sum in Eq. 60 gives $2L^{2}-2L$ independent constraints. The second sum gives $L$ independent constraints. The third sum gives no independent constraints, since all terms in the third sum can be written as a product of terms in the first and second sums. The fourth sum gives $L-1$ independent constraints. The fifth term gives a single independent constraint. The final term gives no independent constraints, since it is a product of the terms in the fourth sum and the fifth term. So, in total we have $2L^{2}$ independent constraints. The ground state is thereby unique.

Appendix B Majorana Mean Field Theory

Here will derive the Majorana mean field theory for the spin model Eq. 11 on a square $L\times L$ lattice. This will be done by decomposing each spin into 4 Majorana fermions, $\gamma^{1}$, $\gamma^{2}$, $\gamma^{3}$, and $\gamma^{4}$. These Majorana fermions obey the normal Majorana algebra, $\{\gamma^{i},\gamma^{j}\}=0$ and $\gamma^{i}\gamma^{i}=1$. In terms of the Majorana fermions, the spin degrees of freedom can be rewritten as $\sigma^{x}=i\gamma^{1}\gamma^{2}$, $\sigma^{y}=i\gamma^{1}\gamma^{3}$, $\sigma^{z}=i\gamma^{1}\gamma^{4}$, with the local constraint that $\gamma^{1}\gamma^{2}\gamma^{3}\gamma^{4}=1$. It is straightforward to verify that the spin operators defined this way anti-commute.

Let us now consider the terms in Eq. 11. First there are the terms proportional to $K$, involving $\sigma^{z}_{\textbf{r}}\sigma^{z}_{\textbf{r}+\hat{y}}$. Due to the aforementioned constraint, $\sigma^{z}=i\gamma^{1}\gamma^{4}=-i\gamma^{2}\gamma^{3}$. We can thereby write the $K$ terms of Eq. 11 as

[TABLE]

Second, there are the terms proportional to $J$ $\prod_{r\in l(x_{o})}\sigma^{x}_{\textbf{r}}\sigma^{x}_{\textbf{r}+\hat{x}}$. In terms of the Majorana fermions $\sigma^{x}=i\gamma^{1}\gamma^{2}=-i\gamma^{3}\gamma^{4}$. We can thereby write the $J$ terms of Eq. 11 as

[TABLE]

where the product is over $l(x_{o})=\{\textbf{r}=(x,y)|x=x_{o},$ and $1\leq y\leq L\}$.

We will now introduce our mean field ansatz. First, note that $\gamma^{3}_{\textbf{r}}\gamma^{4}_{\textbf{r}+\hat{y}}$ and $\gamma^{1}_{\textbf{r}}\gamma^{2}_{\textbf{r}+\hat{y}}$ commute with all terms in the Hamiltonian. The $K$ terms can then be rewritten as

[TABLE]

where $K_{1,2,\textbf{r}+\hat{y},\textbf{r}}\equiv K\langle\gamma^{3}_{\textbf{r}}\gamma^{4}_{\textbf{r}+\hat{y}}\rangle$ and $K_{3,4,\textbf{r},\textbf{r}+\hat{y}}\equiv K\langle\gamma^{1}_{\textbf{r}+\hat{y}}\gamma^{2}_{\textbf{r}}\rangle$. Eq. 63 is minimized by having $\langle\gamma^{3}_{\textbf{r}}\gamma^{4}_{\textbf{r}+\hat{y}}\rangle=\langle\gamma^{1}_{\textbf{r}+\hat{y}}\gamma^{2}_{\textbf{r}}\rangle=1$. This will correspond to a state with $\langle\sigma^{z}_{\textbf{r}}\sigma^{z}_{\textbf{r}+\hat{y}}\rangle=1$ as expected from the spin-model analysis.

More care must be taken with the $J$ terms due to the boundary conditions. If will be useful to rewrite the $J$ terms as

[TABLE]

where the product is over $l^{\prime}(x_{o})=\{\textbf{r}=(x,y)|x=x_{o}$ and $1\leq y<L\}$. We can now rewrite $J$ terms in the mean field form

[TABLE]

In the above equation

[TABLE]

If we redefine $J_{1,2,\textbf{r}+\hat{y},\textbf{r}}\equiv J_{1,2,\textbf{r}+\hat{y},\textbf{r}}+K_{1,2,\textbf{r}+\hat{y},\textbf{r}}$ and $J_{3,4,\textbf{r},\textbf{r}+\hat{y}}\equiv J_{3,4,\textbf{r},\textbf{r}+\hat{y}}+K_{3,4,\textbf{r},\textbf{r}+\hat{y}}$, we can write the full mean field Hamiltonian as

[TABLE]

where $J_{i,j,\textbf{r},\textbf{r'}}$ is non-zero for the following combinations:

$i=1$, $j=2$, $\textbf{r}=(x,y)$ with $0<y\leq L$ and $\textbf{r}^{\prime}=\textbf{r}-\hat{y}$.

$i=3$, $j=4$, $\textbf{r}=(x,y)$ with $0\leq y<L$ and $\textbf{r}^{\prime}=\textbf{r}+\hat{y}$.

$i=4$, $j=1$, $\textbf{r}=(x,0)$ with $0\leq x<L$ and $\textbf{r}^{\prime}=\textbf{r}+\hat{x}$.

$i=3$, $j=2$, $\textbf{r}=(x,L)$ with $0\leq x<L$ and $\textbf{r}^{\prime}=\textbf{r}+\hat{x}$.

Eq. 18 gives the Majorana dimerization pattern shown in Fig. 5. Each term in the mean field Hamiltonian Eq. 18 commute. The ground state will then minimize each of the Majorana bi-linears in Eq. 18. Using the aforementioned identifications, it is clear that states which minimize the Majorana mean field Hamiltonian will also minimize each of the commuting elements in Eq. 11 as desired.

Appendix C Gauging the Subsystem Symmetry

Here, we will consider gauging the subsystem symmetry generated by $G[l(x_{o})]$ in Eq. 12. The paradigm of gauging a symmetry is to make the symmetry local. The local version of the symmetry generated by $G[l(x_{o})]$ is generated by $\sigma^{x}_{\textbf{r}}\sigma^{x}_{{\textbf{r}}+x}$. This term does not commute with the Hamiltonian Eq. 11. In order to make this local transformation a symmetry of the model, we will proceed in the standard fashion of adding in additional degrees of freedom, i.e., adding gauge fields.

We will add additional spin-1/2 degrees of freedom $\tau$ at the center of each plaquette of the lattice as in Fig. 7a. After this we will change the first term of the Hamiltonian to

[TABLE]

where $\textbf{r}\pm x\pm y$ are the $\tau$ spins indicated in Fig. 7b. We can then flip $\sigma^{z}_{\textbf{r}}$ and $\sigma^{z}_{{\textbf{r}}+x}$ provided we also flip the two $\tau$ spins $\tau^{z}_{\textbf{r}+x+y}$ and $\tau^{z}_{\textbf{r}+x-y}$. So in Eq. 68 we have gauged the subsystem symmetry as desired. As the symmetry is now local, there is no need to include the non-local terms proportional to $J$. We can now include a term to energetically enforce invariance under the new local symmetry. This is done via the term

[TABLE]

Let us now redefine the lattice as in Fig. 7c, and relabel $\tau\rightarrow\sigma$ (which should not cause confusion since the $\tau$ and $\sigma$ operators are defined on different lattice sites). After this relabeling we find that this Hamiltonian is

[TABLE]

where the sum is over the elementary plaquettes $p$ and vertices $v$. This is exactly the Hamiltonian for the deconfined $\mathbb{Z}_{2}$ lattice gauge theory, i.e., the toric code. Our spin model can thereby be identified as a $\mathbb{Z}_{2}$ lattice gauge theory where we have ”un-gauged” the $\mathbb{Z}_{2}$ gauge symmetry into a subsystem symmetry.