# Multipole and fracton topological order
via gauging foliated SPT phases

Hiromi Ebisu1 1Yukawa Institute for Theoretical Physics, Kyoto University
Masazumi Honda2 1Yukawa Institute for Theoretical Physics, Kyoto University
Taiichi Nakanishi2,1 1Yukawa Institute for Theoretical Physics, Kyoto
University

###### Abstract

Spurred by recent development of fracton topological phases, unusual
topological phases possessing fractionalized quasi-particles with mobility
constraints, the concept of symmetries has been renewed. In particular, in
accordance with the progress of multipole symmetries, associated with
conservation of multipoles, such as dipole or quadruple moments as well as
global charges, there have been proposed topological phases with such
symmetries. These topological phases are unconventional as excitations are
subject to mobility constraints corresponding to the multipole symmetries. We
demonstrate a way to construct such phases by preparing layers of symmetry
protected topological (SPT) phases and implementing gauging a global symmetry.
After gauging, the statistics of a fractional excitation is altered when
crossing the SPT phases, resulting in topological phases with the multipole
symmetries. The way we construct the phases allows us to have a comprehensive
understanding of field theories of topological phases with the multipole
symmetries and other fracton models.

###### Contents

  1. 1 Introduction
  2. 2 BF theory with dipole symmetry
  3. 3 Gauging SPT phases
    1. 3.1 Dipole symmetry in $x$-direction
    2. 3.2 Alternative derivation of the GSD – Fusion rules
  4. 4 Dipole symmetries in $x$\- and $y$-direction
  5. 5 Construction of other fracton models
    1. 5.1 Exotic $\mathbb{Z}_{N}$ theory
    2. 5.2 X-cube model
  6. 6 Discussion and conclusion
    1. 6.1 Comments on entanglement entropy
    2. 6.2 Future directions

##  1 Introduction

Topologically ordered phases have been one of the central subjects in
condensed matter physics, and provide us a paradigm shift in understating
phases of matter [1, 2, 3, 4]. The most prominent feature of these phases is
that they support fractionalized quasi-particle excitations, the so-called
anyons [5, 6]. Such exotic excitations are not only theoretically and
experimentally intriguing but also are of practically importance for utilizing
quantum computers [7, 8, 9, 10].

Recently, there have been proposed new topological phases, referred to as the
fracton topological phases [11, 12, 13]. The distinct property compared with
the conventional topologically ordered phases is that the fractional
excitations are subject to mobility constraints. Due to this feature, ground
state degeneracy (GSD) becomes sub-extensive. Compared with the conventional
topologically ordered phases, where the GSD depends on topology of a manifold,
the sub-extensive GSD found in the fracton topological phases implies that one
cannot separate the UV and IR physics, which necessitates a new theoretical
framework.

To aim for establishing a consistent continuum field theoretical description
of the fracton topological phases, recently the concept of symmetries has been
reconsidered. One of such update is multipole symmetry, which is a
generalization of ordinary global symmetry in the sense that multipole
moments, such as dipole or quadrupole moments in addition to global charges
are conserved [14, 15, 16, 17, 18, 19]. More explicitly, multipole symmetries
correspond to a case when a theory is invariant under a global phase rotation
which has a polynomial form of the spatial coordinate. For instance, in the
case of a scalar theory with a field $\Phi$, we say the theory respects dipole
symmetry when the theory is invariant under global phase shift in the form
$\Phi\to e^{i(\alpha+x\cdot\beta)}\Phi$, where $\alpha$, $\beta$, and $x$
represent a constant number, a constant spatial vector, and the linear spatial
coordinate vector, respectively. Other cases of higher order multipole
symmetries, such as quadrupole symmetry can be analogously defined, by
considering the phase shift which depends on higher order spatial coordinates.
Motivated by this progress, several topological phases with such a symmetry
were studied [18, 20, 21, 22, 23]. These phases have fractional excitations
which are subject to a mobility constraint, corresponding to the multipole
symmetry. Due to this feature, the phases exhibit unusual GSD dependence on
the system size; the GSD depends on the greatest common divisor between
integer $N$ characterizing the fractional charge and the system size. Although
several gauge theories and topological models with such symmetries have been
recently studied, it is still largely unknown that what physical implication
or what kind of physical properties one can learn when a theory respects such
symmetries.

In this paper, we demonstrate a way to establish topological phases with the
multipole symmetries. In doing so, we prepare stack of symmetry protected
topological (SPT) phases [24, 25, 26] with global $\mathbb{Z}_{N}$ symmetry,
which are invertible gapped phases respecting the global symmetry
$\mathbb{Z}_{N}$, and implement gauging. Generically, it is known that via
gauging, a $d+1$-dimensional SPT phase with global symmetry $G$, characterized
by a $(d+1)$-cohomology class is mapped to the so-called Dijkgraaf-Witten
theory [27] with a twist term, corresponding to the non-trivial cohomology
class [28, 29, 30, 31]. Such a twist term, that we dub the Dijkgraaf-Witten
(DW) twist term throughout this paper, modifies braiding statistics of anyons
in the theory. One example of the DW twist terms defined in $(2+1)$d111
Throughout this paper, we use two letters “d” and “D” to abbreviate the
dimensions. The small “d” stands for space-time dimensions whereas the capital
one “D” does for spatial dimensions.  has the form $\sim a^{1}\wedge
a^{2}\wedge a^{3}$ where $a^{j}\;(j=1,2,3)$ denotes a $1$-form gauge field
associated with a global symmetry [28, 29, 30, 31]. Key insight of our
construction is that via gauging the global symmetry, the stack of the SPT
phases become arrays of the DW twist terms to change statistics of quasi-
particles, which results in a topological phase with the multipole symmetries.
By this approach, one can make intuitive understanding of the BF theories with
multipole symmetries that was previously discussed [32]. To our knowledge, the
present paper is the first work to address how the DW twist term plays a
pivotal role in the context of the fracton topological phases, especially,
topological phases with multipole symmetries.

Our proposal may share the same spirit as an attempt to construct fracton
topological models, such as the X-cube model via defect condensations [33].
Compared with this attempt, our construction is simpler as we do not impose
energetic projection on the models to implement the condensation. Our approach
allows us to make comprehensive understanding of foliated BF theories [34, 35,
36, 37, 38, 32], which are special topological field theories describing
fracton and multipole topological order by interpreting coupling terms
$A^{I}\wedge b\wedge e^{I}$ (to be defined in the next section) as the
foliated DW twist terms. Moreover, our result would comply with growing
interests in unified construction of topological defect lines obtained by
stacking SPT phases and gauging [39], implying a possibility to explore the
multipole topological phases in a broader context. As a byproduct, our
construction of the stacking SPT phases and gauging leads us to the other
fracton model, such as the exotic $\mathbb{Z}_{N}$ gauge theory [40] and
X-cube model [13]. In this regard, the way we construct phases may provide us
a unified insight on topological phases with multipole symmetries and other
fracton models.

The rest of this work is organized as follows. In Sec. 2, to extract an
intuition behind our construction, we introduce a BF theory description of
topological phases with multipole symmetries, following the procedure
discussed in [32]. In Sec. 3, we study a UV lattice model to construct the
topological phase with the multipole symmetry, especially dipole symmetry in
the $x$-direction, which is the simplest case of the multipole symmetry. We
explicitly demonstrate how the DW twist term affects the statistics of
fractional excitations, giving the desired topological phase with the dipole
symmetry. We also generalize the discussion to the other case of the multipole
symmetry in Sec. 4. In Sec. 5, we apply our approach to other fracton models,
the so-called exotic $\mathbb{Z}_{N}$ gauge theory [40] and the $X$-cube model
[13]. Finally, we conclude our work in Sec. 6 with a few remarks and future
research directions.

##  2 BF theory with dipole symmetry

To gain intuition behind our idea to construct topological phases with
multipole symmetry via gauging SPT phases, we begin by reviewing how
topological BF theories with multipole symmetries are introduced, starting
with an argument of these symmetries and gauge fields associated with them.
Throughout this section, we focus on $(2+1)$d system and for the sake of the
notational brevity, we employ the differential forms.

Suppose we have a theory with conserved charges associated with global $U(1)$
and dipole symmetry in the $x$-direction, described by $Q$ and $Q_{x}$
respectively. These charges are subject to the following relation [41, 18, 32]

$[iP_{I},Q]=0\quad[P_{x},Q_{x}]=Q,\quad[P_{y},Q_{x}]=0\qquad(I=x,y).$ (1)

Here, $P_{I}$ denotes the translation operator along the $I$-th direction.
While the first relation is obvious as $Q$ represents the global charge, the
intuition behind the second and third relations is that one associates the
charge $Q_{x}$ with the dipole moment in the $x$-direction, $x\rho$. Indeed,
by writing the charge and dipole moment as $\rho$ and $x\rho$, where $\rho$
denotes the density of the $U(1)$ charge, we think of shifting them by a
constant in $x$\- or $y$-direction. For instance, if one shifts the dipole
moment $x\rho$ by a constant in the $x$-direction, the change of dipole moment
under the shift is given by $(x+\Delta x)\rho-x\rho=\Delta x\rho$, where
$\Delta x$ is constant, corresponding to the nontrivial commutation relation
as described in the second relation in (1). The third relation is analogously
discussed.

Introducing $1$-form currents $j$ and $K_{x}$, we write the charges $Q$ and
$Q_{x}$, via the integral expression using the conserved current as

$Q=\int_{V}*j,\quad Q_{x}=\int_{V}*K_{x}.$

Here, $*$ represents the Hodge dual and $V$ does spatial volume. To reproduce
the relation (1), we demand that

$*K_{x}=*k_{x}-x*j$ (2)

with $k_{x}$ being a local (non-conserved) current. A simple calculation
verifies the relation (1). We also introduce $U(1)$ $1$-form gauge fields $a$,
$A$ with the coupling term defined by

$S_{\rm dipole}=\int_{V}\left(a\wedge*j+A\wedge*k_{x}\right).$ (3)

We need to have an appropriate gauge transformation such that the condition of
the coupling term being gauge invariant yields the conservation law of the
currents. The following gauge transformation does the job:

$a\to a+d\Lambda+\Sigma dx,\quad A\to A+d\Sigma.$ (4)

Here, $\Lambda$ and $\Sigma$ denote the gauge transformation parameters. One
can check that the gauge invariance of the coupling term $S_{\rm dipole}$
under the gauge transformation (4) gives rise to

$d*j=0,\quad d\left(*k_{x}-x*j\right)=d*K_{x}=0.$

Now we introduce a BF theory using the gauge fields of the multipole symmetry.
We define the following gauge invariant fluxes as

$f\vcentcolon=da-A\wedge dx,\quad F\vcentcolon=dA,$ (5)

and put these fluxes in the BF theory format:

$\mathcal{L}_{\rm dipole}=\frac{N}{2\pi}b\wedge f+\frac{N}{2\pi}c\wedge F,$
(6)

with $b$ and $c$ representing $U(1)$ $1$-form gauge fields. We also define the
foliation field, which is introduced in the context of field theory of fracton
topological phases, following [34, 35]. Generally, foliation is defined by co-
dimension one sub-manifolds stacked along a direction, characterized by the
$1$-form field, which we call the foliation field. In the present case, we set
the foliation field by $e^{x}\vcentcolon=dx$. Rewriting (6), we arrive at the
following BF theory:

$\boxed{\mathcal{L}_{\rm dipole}=\frac{N}{2\pi}a\wedge
db+\frac{N}{2\pi}A\wedge dc+\frac{N}{2\pi}A\wedge b\wedge e^{x}.}$ (7)

The BF theory consists of two layers of the toric codes and the couplings
between the two. Generally, the BF theory of a topological phase with
multipole symmetries has the similar form as (7); the BF theory consists of
the layers of the toric codes with couplings between the layers, where the
number of layers corresponds to the number of degrees of freedom of the
multipoles (meaning, the number of global charge, dipole, quadrupole, e.t.c.)
[32]. In addition to the gauge symmetry (4), the BF theory (7) also admits the
following gauge symmetry with respect to $b$ and $c$:

$b\to b+d\lambda,\quad c\to c+d\gamma-\lambda e^{x}.$ (8)

with $\lambda$ and $\gamma$ being gauge parameters.

To extract physical intuitions behind this BF theory, let us rewrite the
theory in terms of tensor gauge fields. One of the approaches for that is to
integrate out the coupling term, corresponding to the last term in (7), which
put a several constraints on the gauge fields. Integrating out $b_{0}$ gives
the following condition:

$\partial_{x}a_{y}-\partial_{y}a_{x}=-A_{y}.$ (9)

We also integrate out the $A_{0}$, and obtain the following condition

$\partial_{x}c_{y}-\partial_{y}c_{x}=b_{y}.$ (10)

One can eliminate the gauge fields $b_{y}$ and $A_{y}$ by substituting the
relations (9)(10) into (7). Furthermore, we introduce gauge fields by

$A_{(xx)}\vcentcolon=\partial_{x}a_{x}-A_{x},\quad
c_{(xx)}\vcentcolon=\partial_{x}c_{x}+b_{x},$ (11)

whose gauge transformation reads [referring to (4)(8)]

$\displaystyle A_{(xx)}\to A_{(xx)}+\partial_{x}^{2}\Lambda,\quad c_{(xx)}\to
c_{(xx)}+\partial_{x}^{2}\gamma.$ (12)

After the substitution, the Lagrangian (7) has the following form

$\displaystyle\mathcal{L}_{\rm dipole}$ $\displaystyle=$
$\displaystyle\frac{N}{2\pi}\biggl{[}a_{0}(\partial_{x}^{2}c_{y}-\partial_{y}c_{(xx)})+A_{(xx)}(\partial_{t}c_{y}-\partial_{y}c_{0})+a_{y}(\partial_{t}c_{(xx)}-\partial_{x}^{2}c_{0})\biggr{]}.$
(13)

The BF theory (13) is reminiscent of the one of the toric code [7] with a
crucial difference being that the spatial derivative operator in the
$x$-direction is replaced with the second order one.

One of the challenges to study the BF theory such as the one in (13) and those
of other fracton models is the presence of the higher order spatial
derivatives, yielding the UV/IR mixing. To circumvent this problem, we follow
an approach alluded in [42, 42, 18]. Instead of directly investigating the
theory (13), we consider the following BF theory defined on a discrete
Euclidean lattice where each gauge field takes integer value (see [42, 42, 18,
32] for more detailed discussion on this point )

$\displaystyle\hat{\mathcal{L}}_{\rm dipole}$ $\displaystyle=$
$\displaystyle\frac{2\pi}{N}\biggl{[}\hat{a}_{0}(\Delta_{x}^{2}\hat{c}_{y}-\Delta_{y}\hat{c}_{(xx)})+\hat{A}_{(xx)}(\Delta_{\tau}\hat{c}_{y}-\Delta_{y}\hat{c}_{0})+\hat{a}_{y}(\Delta_{\tau}\hat{c}_{(xx)}-\Delta_{x}^{2}\hat{c}_{0})\biggr{]}.$
(14)

Here, the gauge fields take values in $\mathbb{Z}_{N}$ and $\Delta_{\tau,x,y}$
denotes the discretized differential operator. This BF theory admits the
following gauge symmetry:

$\displaystyle\hat{a}_{0}\to\hat{a}_{0}+\Delta_{\tau}\lambda_{a},\quad\hat{A}_{(xx)}\to\hat{A}_{(xx)}+\Delta_{x}^{2}\lambda_{a},\quad\hat{a}_{y}\to\hat{a}_{y}+\Delta_{y}\lambda_{a}$
$\displaystyle\hat{c}_{0}\to\hat{c}_{0}+\Delta_{\tau}\lambda_{c},\quad\hat{c}_{(xx)}\to\hat{c}_{(xx)}+\Delta_{x}^{2}\lambda_{c},\quad\hat{c}_{y}\to\hat{c}_{y}+\Delta_{y}\lambda_{c},$
(15)

where $\lambda_{a}$ and $\lambda_{c}$ are integer gauge parameters. Equation
of motions of the theory (14) implies that the following gauge invariant field
strength vanish

$\displaystyle
E_{ax}=\Delta_{\tau}\hat{A}_{(xx)}-\Delta_{x}^{2}\hat{a}_{0},\quad
E_{ay}=\Delta_{\tau}\hat{a}_{y}-\Delta_{y}\hat{a}_{0},\quad
B_{a}=\Delta_{x}^{2}\hat{a}_{y}-\Delta_{y}\hat{A}_{(xx)},$ $\displaystyle
E_{cx}=\Delta_{\tau}\hat{c}_{(xx)}-\Delta_{x}^{2}\hat{c}_{0},\quad
E_{cy}=\Delta_{\tau}\hat{c}_{y}-\Delta_{y}\hat{c}_{0},\quad
B_{c}=\Delta_{x}^{2}\hat{c}_{y}-\Delta_{y}\hat{c}_{(xx)}.$ (16)

The equation of motions ensures that there is no non-trivial local gauge
invariant operators. Yet, the theory (14) does have non-local gauge-invariant
operators, which can be constructed from the gauge fields either
($\hat{a}_{0},\hat{A}_{(xx)},\hat{a}_{y}$) or
($\hat{c}_{0},\hat{c}_{(xx)},\hat{c}_{y}$). Especially, when the theory (14)
is placed on a torus, it admits non-contractible Wilson loops, contributing to
the GSD. To evaluate this number, we set the theory (14) on a torus geometry
by imposing the periodic boundary condition via $\hat{x}\sim\hat{x}+L_{x}$
$\hat{y}\sim\hat{y}+L_{y}$,222 The coordinate $(\hat{x},\hat{y})$ takes
integer number in the unit of the lattice spacing.  i.e., the torus with the
length in the $x/y$-direction being $L_{x}/L_{y}$. Focusing on the Wilson
loops of the gauge fields $(\hat{A}_{(xx)},\hat{a}_{y})$, the distinct number
of such loops which amounts to be the GSD. In the $x$-direction, we have the
following two types of the Wilson loops [22, 32]

$W_{0}(\hat{y})=\exp\biggl{[}\frac{2\pi
i}{N}\sum_{\hat{x}=1}^{L_{x}}\hat{A}_{(xx)}(\hat{x},\hat{y})\biggr{]},\quad
W_{\rm dipole}(\hat{y})=\exp\biggl{[}\frac{2\pi
i}{N}\alpha_{x}\sum_{\hat{x}=1}^{L_{x}}\hat{x}\hat{A}_{(xx)}(\hat{x},\hat{y})\biggr{]},$
(17)

where $\alpha_{x}=\frac{N}{\gcd(N,L_{x})}$ and $\gcd$ stands for the greatest
common divisor. While the first term in (17) corresponds to the standard
Wilson loop which is also found in other conventional topologically ordered
phases, the second term has the unconventional form. It is interpreted as the
“dipole of the Wilson loop" in the sense that the argument of the exponential
is a linear function of the coordinate $\hat{x}$, reminiscent of the
integration of the dipole moment, $x\rho$. This mirrors the fact that the BF
theory respects the dipole symmetry in the $x$-direction (1). From equation of
motions (16), it follows that the two types of loops are deformable in the
$y$-direction, implying that $W_{0}(\hat{y})$ and $W_{\rm dipole}(\hat{y})$ do
not depend on $\hat{y}$. Hence, the non-contractible Wilson loops of the gauge
field $\hat{A}_{(xx)}$ in the $x$-direction are labeled by
$\mathbb{Z}_{N}\times\mathbb{Z}_{\gcd(N,L_{x})}$.

One can evaluate the non-contractible Wilson loops in the $y$-direction, which
are also labeled by the same quantum number,
$\mathbb{Z}_{N}\times\mathbb{Z}_{\gcd(N,L_{x})}$ [22]. Overall, the GSD, which
is equivalent to the number of distinct non-contractible Wilson loops, is
given by

${\rm GSD}=\left[N\times\gcd(N,L_{x})\right]^{2}.$ (18)

Compared with the other fracton models, such as the $X$-cube model [13], which
exhibit the sub-extensive GSD, the theory shows unusual GSD dependence on the
system size due to the dipole symmetry.

##  3 Gauging SPT phases

Now we come to the main part of this paper. In the previous section, based on
the argument on the global and dipole symmetry, we introduced the topological
field theory (7), which respects such symmetries. To further study the
properties of the theory, we integrate out the coupling term, $A\wedge b\wedge
e^{x}$, corresponding to the third term in (7) to simplify the theory to the
one with higher rank derivative terms (13). Instead of integrating out the
coupling term, in this section, we give more physical comprehension on this
term and based on this understanding, we construct a UV lattice model
corresponding to the topological field theory with the dipole symmetry.

The form $A\wedge b$ that enters (7) reminds us of the so-called the DW twist
term [27, 28, 29, 30, 31]. Recalling that the $e^{x}$ is the foliated field,
along which $(1+1)$d sub-manifolds labeled by $(t,y)$ are stacked, it is
tempting to interpret the BF theory (7) as the one described by the two copies
of the standard BF theories with “foliated DW twists", namely, DW twist terms
$A\wedge b$ stacked in the $x$-direction. In the following, we demonstrate a
concrete lattice model to realize such twist terms.

###  3.1 Dipole symmetry in $x$-direction

To this end, we introduce two copies of 2D square lattices. We define integer
coordinate of a node of one lattice (black dot in Fig. 1(a)) and the one of
the other lattice (blue square in Fig. 1(a)) by $(\hat{x},\hat{y})$ in the
unit of the lattice spacing. For illustration purposes, we intentionally shift
a bit the one of the lattice with respect to the other one, as depicted in
Fig. 1(a). Also, we accommodate $\mathbb{Z}_{N}$ qubit at each site of the two
square lattices, as portrayed in Fig. 1(a). We label the Hilbert space and the
$\mathbb{Z}_{N}$ Pauli operator at each site as $\ket{a}_{i}$
($a\in\mathbb{Z}_{N}$) and $\tau^{X/Z}_{i}$, where the subscript $i=1,2$
distinguishes which square lattice the site belongs to. The Pauli operators
act on the Hilbert spaces as

$\tau^{X}_{i}\ket{a}_{i}=\ket{a+1}_{i},\quad\tau^{Z}_{i}\ket{a}_{i}=\omega^{a}\ket{a}_{i}\quad(\omega:=e^{2\pi
i/N}).$ (19)

(a)

(b)

(c)

(d)

Figure 1: (a) Two square lattices (solid black and dashed blue lines). The
links are oriented according to the arrows. For illustrative purposes, we
intentionally shift a bit one of the lattice with respect to the other. (b)
Along a line at $\hat{x}=\hat{x}_{0}$, we implement the CZ gates (21) on the
$\mathbb{Z}_{N}$ qubit states according to (23). Examples of the terms in (23)
are shown by red and blue colors. (c) After gauging, the resulting Hamiltonian
becomes two copies of the $\mathbb{Z}_{N}$ toric codes. Away from the line at
$\hat{x}=\hat{x}_{0}$, the Hamiltonian composes of the vertex and the
plaquette terms, shown in the four panels at the bottom, whereas at
$\hat{x}=\hat{x}_{0}$, the additional Pauli operator is attached to the term
depicted by the two panels on the right. (d) Trajectory of the magnetic charge
$m_{1}$ when crossing the line at $\hat{x}=\hat{x}_{0}$ in the model (30).

In this setting, we define the following Hamiltonian of a paramagnet:

$H_{0}=-\sum_{\mathbf{r}}\Bigl{[}\tau^{X}_{1,\mathbf{r}}+\tau^{X}_{2,\mathbf{r}}\Bigr{]}.$
(20)

Here, we have defined vectors to abbreviate the coordinates by
$\mathbf{r}\vcentcolon=(\hat{x},\hat{y})$. We introduce non-trivial SPT phases
in this model. For this purpose, we start with introducing one $(1+1)d$ SPT
along a $y$-direction at $\hat{x}=\hat{x}_{0}$. We define the generalized
controlled-Z gate $CZ$ and $CZ^{\dagger}$ acting on two $\mathbb{Z}_{N}$ qubit
states. Its action on the states is given by

$CZ\ket{a}\ket{b}=\omega^{ab}\ket{a}\ket{b},\quad
CZ^{\dagger}\ket{a}\ket{b}=\omega^{-ab}\ket{a}\ket{b}\qquad(a,b\in\mathbb{Z}_{N}).$
(21)

At $\hat{x}=\hat{x}_{0}$, we act the CZ gate on two adjacent $\mathbb{Z}_{N}$
qubit states, composed of one state from a site on a square lattice and the
one on the other lattice. More explicitly, we consider the following operator

$U_{CZ}\vcentcolon=\prod_{\begin{subarray}{c}\mathbf{r}\\\
\hat{x}=\hat{x}_{0}\end{subarray}}\Bigl{[}CZ^{\dagger}_{(1,\mathbf{r}),(2,\mathbf{r})}CZ_{(1,\mathbf{r}),(2,\mathbf{r}-{\mathbf{j}_{y}})}\Bigr{]}$
(22)

with $\mathbf{j}_{y}$ being a unit vector in the $y$-direction,
$\mathbf{j}_{y}\vcentcolon=(0,1)$ [for latter use, we also define
$\mathbf{j}_{x}\vcentcolon=(1,0)$], and act it on the Hamiltonian (20).
Pictorially, the CZ gates act on the states along the zigzag path portrayed in
Fig. 1(b). After the operation, Hamiltonian (20) is transformed as

$\displaystyle H_{0}^{\prime}$ $\displaystyle:=$ $\displaystyle
U_{CZ}^{\dagger}H_{0}U_{CZ}$ (23) $\displaystyle=$
$\displaystyle-\sum_{\begin{subarray}{c}\mathbf{r}\\\
\hat{x}\neq\hat{x}_{0}\end{subarray}}\Bigl{[}\tau^{X}_{1,\mathbf{r}}+\tau^{X}_{2,\mathbf{r}}\Bigr{]}-\sum_{\begin{subarray}{c}\hat{y}\\\
\hat{x}=\hat{x}_{0}\end{subarray}}\Bigl{[}\tau^{Z}_{2,\mathbf{r}}\tau^{X}_{1,\mathbf{r}}(\tau^{Z})^{\dagger}_{2,\mathbf{r}-\mathbf{j}_{y}}+(\tau^{Z})^{\dagger}_{1,\mathbf{r}+\mathbf{j}_{y}}\tau^{X}_{2,\mathbf{r}}\tau^{Z}_{1,\mathbf{r}}\Bigr{]}.$

Example of the last two terms, which are generalized cluster states [24, 25,
26], are depicted in Fig. 1(b).

The Hamiltonian (23) respects a $(\mathbb{Z}_{N})^{2}$ global symmetry
generated by the operator

$\prod_{\mathbf{r}}\tau^{X}_{i,\mathbf{r}}\quad(i=1,2),$ (24)

which allows us to gauge it. Below, we follow the procedure given in [29] to
implement gauging the global symmetry.333 See also [13] for the similar
discussion in the context of the fracton topological phases.  Via gauging, one
promotes the global symmetry (24) to a local symmetry. To do so, we first
introduce gauge degree of freedom on each link, described by the Hilbert space
$\ket{\varphi}_{i,\mathbf{r}\pm\mathbf{j}_{x}/2},\ket{\varphi}_{i,\mathbf{r}\pm\mathbf{j}_{y}/2}$,
($i=1,2$, $\varphi\in\mathbb{Z}_{N}$), where $i$ distinguishes the link of the
two lattices, with $\mathbb{Z}_{N}$ Pauli operators
$X_{i,\mathbf{r}\pm\mathbf{j}_{x}/2},Z_{i,\mathbf{r}\pm\mathbf{j}_{x}/2}$,
$X_{i,\mathbf{r}\pm\mathbf{j}_{y}/2},Z_{i,\mathbf{r}\pm\mathbf{j}_{y}/2}$ that
act on the links, satisfying the similar relation as (19). The original
$\mathbb{Z}_{N}$ spin d.o.f on each node can be thought of as a matter field
while the newly introduced one located on each link can be regarded as
$\mathbb{Z}_{N}$ gauge field.444 In the terminology of lattice gauge theory,
$Z$ corresponds to link variable and $X$ does to its conjugate, which is a
discrete analog of an exponentiation of electric field.  We impose the
condition that the physical Hilbert space has to have the trivial eigenstate
of the following Gauss law operators:555 Whether we put the dagger $\dagger$
or not on the last four terms in (25) depends on the orientation of the links
in accordance with Fig. 1(a), namely, we put (do not put) the dagger on a
operator on a link which emanates from (terminates at) the node $\mathbf{r}$.

$\displaystyle G_{i,\mathbf{r}}=\tau^{X}_{i,\mathbf{r}}\times
X^{\dagger}_{i,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}X^{\dagger}_{i,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}X_{i,\mathbf{r}-\frac{\mathbf{j}_{x}}{2}}X_{i,\mathbf{r}-\frac{\mathbf{j}_{y}}{2}}\quad(i=1,2),$
(25)

that is, physical states $\ket{\rm phys}$ are subject to
$G_{i,\mathbf{r}}\ket{\rm phys}=\ket{\rm phys}$. The operator (25) corresponds
to the local $\mathbb{Z}_{N}$ spin flip symmetry. We also minimally couple the
quadratic spin coupling terms in the original Hamiltonian (23) to the
$\mathbb{Z}_{N}$ gauge field as

$\displaystyle\tau^{Z}_{2,\mathbf{r}}(\tau^{Z})^{\dagger}_{2,\mathbf{r}-\mathbf{j}_{y}}\to\tau^{Z}_{2,\mathbf{r}}Z^{\dagger}_{2,\mathbf{r}-\frac{\mathbf{j}_{y}}{2}}(\tau^{Z})^{\dagger}_{2,\mathbf{r}-\mathbf{j}_{y}},\quad(\tau^{Z})^{\dagger}_{1,\mathbf{r}+\mathbf{j}_{y}}\tau^{Z}_{1,\mathbf{r}}\to(\tau^{Z})^{\dagger}_{1,\mathbf{r}+\mathbf{j}_{y}}Z_{1,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}\tau^{Z}_{1,\mathbf{r}}$
(26)

so that these terms commute with the Gauss law (25).

In the restricted physical state, satisfying $G_{i,\mathbf{r}}=1$, the
original Hamiltonian (23) now becomes

$\displaystyle\tilde{H}_{0}-=\sum_{\begin{subarray}{c}\mathbf{r}\\\
\hat{x}\neq\hat{x}_{0}\end{subarray}}\Bigl{[}V_{1,\mathbf{r}}+V_{2,\mathbf{r}}\Bigr{]}-\sum_{\begin{subarray}{c}\hat{y}\\\
\hat{x}=\hat{x}_{0}\end{subarray}}\Bigl{[}V_{1,\mathbf{r}}\tau^{Z}_{2,\mathbf{r}}Z^{\dagger}_{2,\mathbf{r}-\frac{\mathbf{j}_{y}}{2}}(\tau^{Z})^{\dagger}_{2,\mathbf{r}-\mathbf{j}_{y}}+V_{2,\mathbf{r}}(\tau^{Z})^{\dagger}_{1,\mathbf{r}+\mathbf{j}_{y}}Z_{1,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}\tau^{X}_{2,\mathbf{r}}\tau^{Z}_{1,\mathbf{r}}\Bigr{]},$
(27)

where

$\displaystyle
V_{i,\mathbf{r}}\vcentcolon=X_{i,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}X_{i,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}X^{\dagger}_{i,\mathbf{r}-\frac{\mathbf{j}_{x}}{2}}X^{\dagger}_{i,\mathbf{r}-\frac{\mathbf{j}_{y}}{2}}.$
(28)

Since $\tau^{Z}_{i,\mathbf{r}}$ commutes with (27), we set
$\tau^{Z}_{i,\mathbf{r}}=1$, which gives rise to Hamiltonian described by only
the gauge fields. We further add the following term to the Hamiltonian

$\displaystyle-\sum_{i,\mathbf{p}}P_{i,\mathbf{p}},\quad
P_{i,\mathbf{p}}\vcentcolon=Z_{i,\mathbf{p}-\frac{\mathbf{j}_{y}}{2}}Z_{i,\mathbf{p}+\frac{\mathbf{j}_{x}}{2}}Z^{\dagger}_{i,\mathbf{p}+\frac{\mathbf{j}_{y}}{2}}Z^{\dagger}_{i,\mathbf{p}-\frac{\mathbf{j}_{x}}{2}},$
(29)

which ensures the fluxless condition to make the gauge theory dynamically
trivial. Also, we have defined the coordinate of a plaquette by
$\mathbf{p}\vcentcolon=(\hat{x}+\frac{1}{2},\hat{y}+\frac{1}{2})$. Overall,
the gauged Hamiltonian reads

$\displaystyle H_{\rm gauged}=-\sum_{\begin{subarray}{c}\mathbf{r}\\\
\hat{x}\neq\hat{x}_{0}\end{subarray}}[V_{1,\mathbf{r}}+V_{2,\mathbf{r}}]-\sum_{\mathbf{p}}[P_{1,\mathbf{p}}+P_{2,\mathbf{p}}]-\sum_{\begin{subarray}{c}\mathbf{r}\\\
\hat{x}=\hat{x}_{0}\end{subarray}}[\widetilde{V}_{1,\mathbf{r}}+\widetilde{V}_{2,\mathbf{r}}]+h.c.,$
(30)

where $V_{i,\mathbf{r}}$ and $P_{i,\mathbf{p}}$ are given in (28)(29) and

$\displaystyle\widetilde{V}_{1,\mathbf{r}}\vcentcolon=Z^{\dagger}_{2,\mathbf{r}-\frac{\mathbf{j}_{y}}{2}}X_{1,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}X_{1,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}X^{\dagger}_{1,\mathbf{r}-\frac{\mathbf{j}_{x}}{2}}X^{\dagger}_{1,\mathbf{r}-\frac{\mathbf{j}_{y}}{2}},$
$\displaystyle\widetilde{V}_{2,\mathbf{r}}\vcentcolon=Z_{1,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}X_{2,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}X_{2,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}X^{\dagger}_{2,\mathbf{r}-\frac{\mathbf{j}_{x}}{2}}X^{\dagger}_{2,\mathbf{r}-\frac{\mathbf{j}_{y}}{2}}.$
(31)

These terms are shown in Fig. 1(c). It is straightforward to check individual
term that constitutes the Hamiltonian (30) commutes with one another. The
gauged Hamiltonian (30) describes two copies of the $\mathbb{Z}_{N}$ toric
code [7], composed of the familiar form of the vertex and plaquette terms,
$V_{i,\mathbf{r}},P_{i,\mathbf{p}}\;(i=1,2)$ except along the 1D line at
$\hat{x}=\hat{x}_{0}$, where the a single $Z_{1}$ or $Z_{2}$ operator is
attached to the vertex term (31).

To see the difference between the standard $\mathbb{Z}_{N}$ toric code and our
gauged Hamiltonian, let us investigate the behavior of the excitations of the
model (30). Since the individual term entering the Hamiltonian commutes with
each other, the ground state $\ket{\Omega}$ satisfies

$\displaystyle
V_{i,\mathbf{r}}\ket{\Omega}=\widetilde{V}_{i,\mathbf{r}}\ket{\Omega}=\ket{\Omega}\quad\forall\mathbf{r},i,$
$\displaystyle{P}_{i,\mathbf{p}}\ket{\Omega}=\ket{\Omega}=\ket{\Omega}\quad\forall\mathbf{p},i.$
(32)

To discuss the behavior of quasi-particle excitations in the model, we denote
a fractional excitation carrying an electric charge which violates the first
condition in (32) with eigenvalue $\omega$ as $e_{i,\mathbf{r}}$. Likewise, we
label a fractional excitation carrying an magnetic charge, violating the
second condition in (32) with eigenvalue being $\omega$ by $m_{i,\mathbf{p}}$.
Further, for latter convenience, we put the symbol $\overline{(\cdots)}$ on
the top of a charge to denote its conjugate. For instance, we write an
electric charge violating the first condition in (32) with eigenvalue
$\omega^{-1}$ as $\overline{e}_{1,\mathbf{r}}$.

Away from $\hat{x}=\hat{x}_{0}$, we think of applying a single operator at a
vertical link $X_{1,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}$
($\hat{x}\ll\hat{x}_{0}$) on the ground state, giving a pair of magnetic
charges. Let us write the magnetic charge schematically as $m_{1}$, omitting
the coordinate dependence for simplicity. One can stretch the magnetic charge
by successive action of the $X_{1}$ operators on the ground state in the
horizontal direction. However, when crossing at $\hat{x}=\hat{x}_{0}$ the
magnetic charge $m_{1}$ is dressed by a electric charge $e_{2}$ as
$X_{1,(\hat{x}_{0},\hat{y})}$ does not commute with
$\widetilde{V}_{2,\mathbf{r}}|_{\hat{x}=\hat{x}_{0}}$, see also Fig 1(d).
Analogous line of thought indicates that the magnetic charge $m_{2}$ induced
by the operator $P_{2,\mathbf{p}}$ is dressed with the electric charge
$\overline{e}_{1}$ when traveling across the line at $\hat{x}=\hat{x}_{0}$.
Also, it is straightforward to see that electric charges are intact by the
terms at $\hat{x}=\hat{x}_{0}$. Therefore, we find that

$m_{1}\Rightarrow m_{1}e_{2},\quad m_{2}\Rightarrow
m_{2}\overline{e}_{1},\quad e_{1}\Rightarrow e_{1},\quad e_{2}\Rightarrow
e_{2}.$ (33)

Here, $\Rightarrow$ represents the fractional excitations crossing the line at
$\hat{x}=\hat{x}_{0}$. This implies that the terms at $\hat{x}=\hat{x}_{0}$
(31) modify the statistics of the magnetic charges (see also [43] for the
relevant discussion in a different context.). More concise discussion of this
point based on fusion rules of fractional charges will be presented in the
next subsection.

(a)

(b)

Figure 2:  (a) The same two lattices as Fig 1(a) with the CZ gates (21) being
acted along the 1D zigzag grey dashed lines. (b) Trajectory of the magnetic
excitation $m_{1}$ dressed with the electric charges, $e_{2}$’s in the model
described by the Hamiltonian (34).

After having seen the effect of the SPT phase, now we prepare arrays of such
$(1+1)d$ SPT phases, which form stacked layers. We place the SPT phases so
that each SPT goes along the $y$-direction and is located adjacent with one
another in the $x$-direction. To this end, starting with the paramagnet
Hamiltonian (20), we think of implementing the CZ gates defined in (21) on the
states along the zigzag path demonstrated in Fig. 2(a),666 Similarly to Fig.
1(a), in the actual configuration of Fig. 2(a), nodes of the two lattices are
located at the same position. However, we intentionally shift a bit one of the
lattice with respect to the other for the sake of visual illustration. and
gauge the global $\mathbb{Z}_{N}$ symmetry. A simple modification of the
previous argument shows that the gauged Hamiltonian becomes

$\displaystyle H_{\rm
dipole}=-\sum_{\mathbf{r}}[\widetilde{V}_{1,\mathbf{r}}+\widetilde{V}_{2,\mathbf{r}}]-\sum_{\mathbf{p}}[P_{1,\mathbf{p}}+P_{2,\mathbf{p}}]+h.c.,$
(34)

where each term is given in (29) and (31).

Based on the discussion around (33), one can see that every time the magnetic
charge $m_{1}$ and $m_{2}$ go in the $x$-direction by one unit lattice
spacing, they are dressed by an electric charge, giving

$\displaystyle
m_{1}\xRightarrow{x}m_{1}e_{2}\xRightarrow{x}m_{1}e_{2}^{2}\xRightarrow{x}m_{1}e_{2}^{3}\xRightarrow{x}\cdots\xRightarrow{x}m_{1}{e}_{2}^{N}=m_{1}\xRightarrow{x}m_{1}e_{2}\xRightarrow{x}\cdots,$
$\displaystyle
m_{2}\xRightarrow{x}m_{2}\overline{e}_{1}\xRightarrow{x}m_{2}\overline{e}_{1}^{2}\xRightarrow{x}m_{2}\overline{e}_{1}^{3}\xRightarrow{x}\cdots\xRightarrow{x}m_{2}\overline{e}_{1}^{N}=m_{2}\xRightarrow{x}m_{2}\overline{e}_{1}\xRightarrow{x}\cdots,$
(35)

where equality in the middle comes from the fact that the electric charge is
$\mathbb{Z}_{N}$ and the arrow $\xRightarrow{x}$ represents the fractional
excitation moving in the $x$-direction by a unit of the lattice spacing. The
trajectory of the magnetic charge is the reminiscent of the form of the dipole
of the Wilson loop discussed in the previous subsection (17) as the intensity
of the electric charges attached with the magnetic charge increases one by one
when propagating in the horizontal direction (see also the configuration in
the bottom of Fig. 2(b)).

To see this point more explicitly and how the lattice model (34) corresponds
to the BF theory with the dipole symmetry (7), we consider the model (34) on a
torus geometry with length being $L_{x}$ and $L_{y}$ in the $x$\- and
$y$-direction to evaluate the GSD. The model (34) admits the following non-
contractible Wilson loops in the $x$-direction: (recall that
$\alpha_{x}=\frac{N}{\gcd(N,L_{x})}$)

$\displaystyle
W_{m_{1}e_{2}}(\hat{y})=\biggl{[}\prod_{\hat{x}=1}^{L_{x}}X_{1,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}Z^{\hat{x}}_{2,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}\biggr{]}^{\alpha_{x}},\quad
W_{m_{2}\overline{e}_{1}}(\hat{y})=\biggl{[}\prod_{\hat{x}=1}^{L_{x}}X_{2,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}\left(Z_{1,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}^{\hat{x}}\right)^{\dagger}\biggr{]}^{\alpha_{x}}.$
(36)

These represent the trajectory of the magnetic charges, dressed with electric
charges whose intensity depends on $\hat{x}$. Physical interpretation on these
loops (36) is that in order for the trajectory of the magnetic charge, which
is dressed with electric charges, to form a non-contractible loop, it has to
be compatible with the periodic boundary condition in the $x$-direction; the
configuration of the charge at $\hat{x}$ and the one at $\hat{x}+L_{x}$ has to
be identical ($\text{mod}\mathbb{Z}_{N}$). Thus, to meet this condition, the
magnetic charge has to wind around the torus multiple times (more explicitly,
$\alpha_{x}$ times) in the $x$-direction.

Since the fractional statistics of the electric charges is not altered by the
arrays of the DW twist terms, we have the following Wilson loops of the
electric charges:

$\displaystyle
W_{e_{1}}(\hat{y})=\prod_{\hat{x}=1}^{L_{x}}X_{1,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}},\quad
W_{e_{2}}(\hat{y})=\prod_{\hat{x}=1}^{L_{x}}X_{2,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}.$
(37)

It is straightforward to check these loops (36) and (37) are deformable in the
$y$-direction, hence $\hat{y}$ independent. Since the loops (36) are labeled
by the quantum number $(\mathbb{Z}_{\gcd(N,L_{x})})^{2}$ and the ones (37) by
$(\mathbb{Z}_{N})^{2}$, the GSD, which amounts to be the number of distinct
non-contractible loops of the electric and magnetic charges in the
$x$-direction, is given by

${\rm GSD}=[N\times\gcd(N,L_{x})]^{2},$ (38)

which is the identical to (18) of the BF theory (7) with dipole symmetry.

To summarize the discussion, motivated by the form of the BF theory of a
topological phase with dipole symmetry (7), we have constructed the
corresponding UV lattice model by deciphering the coupling term $A\wedge
b\wedge e^{x}$ as the foliated DW twist terms. In the next subsection, we will
see how such DW twist terms affect the statistics of quasi-particles by
thoroughly investigating the fusion rules.

###  3.2 Alternative derivation of the GSD – Fusion rules

One can derive the same GSD of the model (38) in a different approach by
counting the superselection sectors of the fractional excitations (i.e., the
number of distinct fractional excitations) and making use of the property of
(35), which will be useful for evaluation of GSD in the other topological
phase with multipole symmetries.

When we act a single $Z_{1}$ operator on the ground state at a horizontal link
or vertical link, i.e.,
$Z_{1,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}\ket{\Omega}$, or
$Z_{1,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}\ket{\Omega}$, it induces a pair of
$\mathbb{Z}_{N}$ electric charges, described by the following fusion rule:

$I\to e_{1,\mathbf{r}}\otimes\overline{e}_{1,\mathbf{r}+\mathbf{j}_{x}},\quad
I\to e_{1,\mathbf{r}}\otimes\overline{e}_{1,\mathbf{r}+\mathbf{j}_{y}},$ (39)

Here, $I$ denotes the trivial sector. From this fusion rule, one can identify
the electric charges which are adjacent with one another, namely,

$e_{1,\mathbf{r}}\simeq e_{1,\mathbf{r}+\mathbf{j}_{x}}\simeq
e_{1,\mathbf{r}+\mathbf{j}_{y}}$ (40)

with $\simeq$ representing the identification of the fractional excitations.
Also, due to the periodic boundary condition, we have

$e_{1,\mathbf{r}+L_{x}\mathbf{j}_{x}}\simeq
e_{1,\mathbf{r}+L_{y}\mathbf{j}_{y}}\simeq e_{1,\mathbf{r}}.$ (41)

The two conditions (40) and (41) imply that the electric charge
$e_{1,\mathbf{r}}$ does not depend on the $x,y$ coordinate, allowing us to
set, e.g., $\mathbf{r}=(1,1)\vcentcolon=\mathbf{r}_{0}$,

$e_{1,\mathbf{r}}\simeq e_{1,\mathbf{r}_{0}}(\vcentcolon=\mathfrak{E}_{1}).$
(42)

Likewise, it can be shown that the electric charge $e_{2,\mathbf{r}}$ does not
depend on the $x,y$ coordinate, and one can set $\mathbf{r}=\mathbf{r}_{0}$,

$e_{2,\mathbf{r}}\simeq e_{2,\mathbf{r}_{0}}(\vcentcolon=\mathfrak{E}_{2}).$
(43)

Therefore, the theory (34) admits $(\mathbb{Z}_{N})^{2}$ electric charges,
labeled by $\mathfrak{E}_{1}^{u}\mathfrak{E}_{2}^{v}\;(u,v\in\mathbb{Z}_{N})$.

When it comes to the magnetic charges, we think of applying a single $X_{1}$
operator on the ground state at a horizontal or vertical link, that is,
$X_{1,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}\ket{\Omega}$, or
$X_{1,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}\ket{\Omega}$. In the former case,
we have the following fusion rule

$I\to m_{1,\mathbf{p}}\otimes\overline{m}_{1,\mathbf{p}-\mathbf{j}_{y}},$ (44)

whereas in the latter case, we find

$I\to{m}_{1,\mathbf{p}-\mathbf{j}_{x}}\otimes\overline{m}_{1,\mathbf{p}}\otimes\overline{e}_{2,\mathbf{r}}.$
(45)

The fusion rule (45) is the concise description of what we have seen in the
previous section, that is, when a magnetic charge goes in the $x$-direction by
a unit of the lattice spacing, an electric charge is attached to it (35). From
the fusion rule (44), we have

$m_{1,\mathbf{p}}\simeq{m}_{1,\mathbf{p}-\mathbf{j}_{y}}.$ (46)

From this relation and the periodic boundary condition in the $y$-direction,
$m_{1,\mathbf{p}+L_{y}\mathbf{j}_{y}}\simeq m_{1,\mathbf{p}}$, it follows that
the magnetic charge $m_{1,\mathbf{p}}$ does not depend on the $y$-coordinate,
allowing us to set $\hat{y}=1$, that is,
$\mathbf{p}=(\hat{x}+\frac{1}{2},\frac{3}{2})\vcentcolon=\mathbf{q}_{0}$.
Also, from the fusion rule (45), jointly with (43), one finds

${m}_{1,\mathbf{q}_{0}}\simeq
m_{1,\mathbf{q}_{0}-\mathbf{j}_{x}}\otimes{\overline{\mathfrak{E}}}_{2}.$ (47)

The adjacent magnetic charges in the $x$-direction are identified up to an
electric charge. Similarly to the argument in the $y$-direction, one naively
wonders the magnetic charge also does not depend on $x$ coordinate. However,
we have to carefully take into account the periodic boundary condition in the
$x$-direction and the effect of the electric charge which is attached to the
magnetic charge. Indeed, by iterative use of (47), we have

$m_{1,\mathbf{q}_{0}+K_{x}\mathbf{j}_{x}}\simeq
m_{1,\mathbf{q}_{0}}\otimes(\overline{{\mathfrak{E}}_{2}})^{K_{x}}\;(K_{x}\in\mathbb{Z}).$
(48)

The condition (48) needs to be compatible with the periodic boundary condition
in the $x$-direction. Hence, depending on $N$ and $L_{x}$, the admissible
magnetic excitation is the multiple magnetic charges rather than a single
magnetic charge. We think of the multipole magnetic charges, denoted by
$m^{\alpha_{x}}_{1,\mathbf{q}_{0}}$, with
$\alpha_{x}=\frac{N}{\gcd(N,L_{x})}$, which satisfies

${m}^{\alpha_{x}}_{1,\mathbf{q}_{0}}\simeq
m^{\alpha_{x}}_{1,\mathbf{q}_{0}-\mathbf{j}_{x}}\otimes\overline{{\mathfrak{E}}}^{\alpha_{x}}_{2}$
(49)

and the periodic boundary condition

$m^{\alpha_{x}}_{1,\mathbf{q}_{0}+L_{x}\mathbf{j}_{x}}\simeq
m^{\alpha_{x}}_{1,\mathbf{q}_{0}}.$ (50)

Hence, from these two conditions (49) and (50), it follows that the magnetic
charge $m^{\alpha_{x}}_{1,\mathbf{q}_{0}}$ at any $x$ coordinate is generated
by an electric charge $\mathfrak{E}_{2}$ and
[$\mathbf{p}_{0}\vcentcolon=(\frac{3}{2},\frac{3}{2})$]

$m^{\alpha_{x}}_{1,\mathbf{p}_{0}}(\vcentcolon=\mathfrak{M}_{1}),$ (51)

indicating that the magnetic charge induced by the Pauli operators
$X_{1}^{\prime}$s, is labeled by the $\mathfrak{M}_{1}$. Analogous analysis
leads us to that the magnetic charge induced by the Pauli operators
$X_{2}^{\prime}$s, is labeled by

$m^{\alpha_{x}}_{2,\mathbf{p}_{0}}(\vcentcolon=\mathfrak{M}_{2}).$ (52)

Therefore, the theory supports $(\mathbb{Z}_{\gcd(N,L_{x})})^{2}$ magnetic
charges, which are characterized by
$(\mathfrak{M}_{1})^{w}(\mathfrak{M}_{2})^{s}$ with
$w,s\in\mathbb{Z}_{\gcd(N,L_{x})}$.

In summary, via a close investigation of the fusion rules, we have evaluated
the number of superselection sectors of the model (34). There are $N^{2}$
electric charges and $\gcd(N,L_{x})^{2}$ magnetic charges. Hence, the GSD on
the torus, which amounts to be the number of the superselection sectors, is
given by ${\rm GSD}=[N\times\gcd(N,L_{x})]^{2}$, which is the same as (38).

##  4 Dipole symmetries in $x$\- and $y$-direction

(a)

(b)

(c)

(d)

(e)

Figure 3: (a) Three 2D square lattices. Two of them are introduced so that
their nodes (red squares and green triangles) are located at the center of the
plaquette of the other lattice (black nodes and links). (b) CZ gates (54) are
implemented along the red and green zigzag path. (c) After gauging the global
$\mathbb{Z}_{N}$ symmetry, gauge fields are located on each link (black dot,
red square, and green triangle). (d) The terms defined on the lattice (c) that
constitutes the Hamiltonian (57) and (58). (e) Schematic illustration how the
magnetic charge $m_{1}$ is dressed with electric charges when it goes in the
$x$\- or $y$-direction.

One can generalize the previous argument given in sec. 3.1 to the case where
we prepare stacking of the SPT phases in the $y$-direction as well as the
$x$-direction. Consider three layers of 2D square lattices, The orientation of
the links is the same as the previous case (Fig. 1(a)). We portray the
lattices in Fig. 3(a) with intentionally shifting two of them with respect to
the other for visual illustration. We define the Hilbert space of the
$\mathbb{Z}_{N}$ qubit on each node of the lattices as $\ket{a}_{i}$
$(a\in\mathbb{Z}_{N})$ with index $i$ distinguishing the three lattices.
Pictorially $i=1,2$ and 3 in Fig. 3(a) correspond to the lattice with black,
red and green colors respectively. Then we introduce the Hamiltonian of the
paramagnet as

$\widetilde{H}_{0}=-\sum_{\mathbf{r}}\left(\tau^{X}_{1,\mathbf{r}}+\tau^{X}_{2,\mathbf{r}}+\tau^{X}_{3,\mathbf{r}}\right).$
(53)

Here, $\tau_{i,\mathbf{r}}^{X}$ is the spin flip operator at a node of the
three lattices. We implement the CZ gates on the lattices along the $x$
(green)- and $y$ (red)-direction in accordance with the zigzag path in Fig.
3(b). More explicitly, we act the following operator consisting of the CZ
gates

$\displaystyle\widetilde{U}_{CZ}\vcentcolon=\prod_{\mathbf{r}}CZ^{\dagger}_{(1,\mathbf{r}),(2,\mathbf{r})}CZ_{(1,\mathbf{r}),(2,-\mathbf{j}_{y})}\times
CZ^{\dagger}_{(1,\mathbf{r}),(3,\mathbf{r})}CZ_{(1,\mathbf{r}),(3,\mathbf{r}-\mathbf{j}_{x})}.$
(54)

on the Hamiltonian (53) which yields

$\displaystyle\widetilde{U}^{\dagger}_{CZ}\widetilde{H}_{0}\widetilde{U}_{CZ}=-\sum_{\mathbf{r}}\biggl{[}(\tau^{Z})^{\dagger}_{3,\mathbf{r}-\mathbf{j}_{x}}(\tau^{Z})^{\dagger}_{2,\mathbf{r}-\mathbf{j}_{y}}\tau^{X}_{1,\mathbf{r}}\tau^{Z}_{2,\mathbf{r}}\tau^{Z}_{3,\mathbf{r}}+(\tau^{Z})^{\dagger}_{1,\mathbf{r}+\mathbf{j}_{y}}\tau^{X}_{2,\mathbf{r}}\tau^{Z}_{1,\mathbf{r}}+Z^{\dagger}_{1,\mathbf{r}+\mathbf{j}_{x}}\tau^{X}_{3,\mathbf{r}}\tau^{Z}_{1,\mathbf{r}}\biggr{]}.$
(55)

The Hamiltonian (55) respects $(\mathbb{Z}_{N})^{3}$ global symmetry,
corresponding to global spin flip on the three layers of the lattices,

$\prod_{\mathbf{r}}\tau^{X}_{i,\mathbf{r}}\quad(i=1,2,3),$

which we can gauge. Following the same procedure explained in the previous
section (sec. 3.1), gauged Hamiltonian is given by

$\widetilde{H}_{\rm
dipole}=-\biggl{[}\sum_{\mathbf{r}}{V}_{1,\mathbf{r}}+\sum_{\mathbf{p}}{P}_{1,\mathbf{p}}\biggr{]}-\biggl{[}\sum_{\mathbf{r}}{V}_{2,\mathbf{r}}+\sum_{\mathbf{p}}{P}_{2,\mathbf{p}}\biggr{]}-\biggl{[}\sum_{\mathbf{r}}{V}_{3,\mathbf{r}}+\sum_{\mathbf{p}}{P}_{3,\mathbf{p}}\biggr{]}+h.c.$
(56)

with

$\displaystyle{V}_{1,\mathbf{r}}$ $\displaystyle\vcentcolon=$ $\displaystyle
Z^{\dagger}_{3,\mathbf{r}-\frac{\mathbf{j}_{x}}{2}}Z^{\dagger}_{2,\mathbf{r}-\frac{\mathbf{j}_{y}}{2}}X_{1,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}X_{1,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}X^{\dagger}_{1,\mathbf{r}-\frac{\mathbf{j}_{x}}{2}}X^{\dagger}_{1,\mathbf{r}-\frac{\mathbf{j}_{y}}{2}},$
$\displaystyle{V}_{2,\mathbf{r}}$ $\displaystyle\vcentcolon=$ $\displaystyle
Z_{1,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}X_{2,\mathbf{r}-\frac{\mathbf{j}_{y}}{2}}^{\dagger}X_{2,\mathbf{r}-\frac{\mathbf{j}_{x}}{2}}^{\dagger}X_{2,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}X_{2,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}},$
$\displaystyle{V}_{3,\mathbf{r}}$ $\displaystyle\vcentcolon=$ $\displaystyle
Z_{1,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}X_{3,\mathbf{r}-\frac{\mathbf{j}_{y}}{2}}^{\dagger}X_{3,\mathbf{r}-\frac{\mathbf{j}_{x}}{2}}^{\dagger}X_{3,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}X_{3,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}},$
(57)

and

$\displaystyle
P_{i,\mathbf{p}}\vcentcolon=Z_{i,\mathbf{p}-\frac{\mathbf{j}_{y}}{2}}Z_{i,\mathbf{p}+\frac{\mathbf{j}_{x}}{2}}Z^{\dagger}_{i,\mathbf{p}+\frac{\mathbf{j}_{y}}{2}}Z^{\dagger}_{i,\mathbf{p}-\frac{\mathbf{j}_{x}}{2}},$
(58)

where $X_{i}/Z_{i}$ represents the Pauli operator of the gauge field which
reside on the links (Fig. 3(c)). These terms are depicted in Fig. 3(d). The
Hamiltonian (56) with (57) and (58) resembles three copies of the
$\mathbb{Z}_{N}$ toric codes with the difference being that a few Pauli
operators are multiplied with the vertex terms (57), which comes from the
foliated DW twist terms. The ground state of the model (56) is the projected
state $\ket{\Omega}$, satisfying

$\displaystyle{V}_{i,\mathbf{r}}\ket{\Omega}=\ket{\Omega},\quad
P_{i,\mathbf{p}}\ket{\Omega}=\ket{\Omega}\qquad\forall
i,\mathbf{r},\mathbf{p}.$ (59)

By the similar argument in Sec. 3.2, one can discuss how behavior of the
fractional excitations is affected by the foliated DW twist terms. To this
end, we denote a fractional electric charge which violates the first condition
in (59) with eigenvalue $\omega$ by $e_{i,\mathbf{r}}$. Similarly, we label a
fractional magnetic charge which violates the second condition given in (59)
with eigenvalue $\omega$ by $m_{i,\mathbf{p}}$.

It can be shown that no additional charge is attached to electric charges when
they move. Regarding the magnetic charges, depending on the direction they go,
electric charges are attached to them. For instance, by the same logic given
in around (35), one can see that the magnetic charge $m_{1}$ (we omit the
coordinate for simplicity) is dressed with an electric charge $e_{2}$
($e_{3}$) every time it travels in the $x(y)$-direction, schematically
described by (see also Fig. 3(e))

$\displaystyle
m_{1}\xRightarrow{x}m_{1}e_{2}\xRightarrow{x}m_{1}e_{2}^{2}\xRightarrow{x}m_{1}e_{2}^{3}\xRightarrow{x}\cdots\xRightarrow{x}m_{1}{e}_{2}^{N}=m_{1}\xRightarrow{x}m_{1}e_{2}\xRightarrow{x}\cdots,$
$\displaystyle
m_{1}\xRightarrow{y}m_{1}{e}_{3}\xRightarrow{y}m_{1}{e}_{3}^{2}\xRightarrow{y}m_{1}{e}_{3}^{3}\xRightarrow{y}\cdots\xRightarrow{y}m_{1}{e}_{3}^{N}=m_{1}\xRightarrow{y}m_{1}{e}_{3}\xRightarrow{y}\cdots.$
(60)

Here $\xRightarrow{x/y}$ indicates the fractional excitation goes in the
$x/y$-direction by a unit of the lattice spacing. By the same token, one finds
that whenever the magnetic charge $m_{2}$ ($m_{3}$) goes in the
$x$($y$)-direction, an electric charge is attached to it, described by

$\displaystyle
m_{2}\xRightarrow{x}m_{2}\overline{e}_{1}\xRightarrow{x}m_{2}\overline{e}_{1}^{2}\xRightarrow{x}m_{2}\overline{e}_{1}^{3}\xRightarrow{x}\cdots\xRightarrow{x}m_{2}\overline{e}_{1}^{N}=m_{2}\xRightarrow{x}m_{2}\overline{e}_{1}\xRightarrow{x}\cdots,$
$\displaystyle
m_{3}\xRightarrow{y}m_{3}\overline{e}_{1}\xRightarrow{y}m_{3}{e}_{1}^{2}\xRightarrow{y}m_{3}\overline{e}_{1}^{3}\xRightarrow{y}\cdots\xRightarrow{y}m_{3}\overline{e}_{1}^{N}=m_{3}\xRightarrow{y}m_{3}\overline{e}_{1}\xRightarrow{y}\cdots.$
(61)

Also, one can verify that no additional charge is attached to the magnetic
charge $m_{2}$ ($m_{3}$) when it travels in the $y$-($x$-)direction.

We can evaluate the GSD of this model on the torus geometry by investigating
the superselection sectors of the fractional excitations, analogously to Sec.
3.2. By considering the fusion rules of the electric charges, obtained by
acting the Pauli operator $Z_{1}$, $Z_{2}$ and $Z_{3}$ on the ground state,
jointly with the periodic boundary condition in the $x$\- and $y$-direction,
one finds that the electric charges, $e_{i,\mathbf{r}}$ at any coordinate of
$\hat{x}$ and $\hat{y}$, are labeled by

$e_{1,\mathbf{r}}\simeq e_{i,\mathbf{r}_{0}}(\vcentcolon=\mathfrak{E}_{i}),$
(62)

where $i=1,2,3$ and $\mathbf{r}_{0}=(1,1)$. Thus, there are $N^{3}$ electric
charges. As for the magnetic charges, by acting the Pauli operator $X_{1}$ on
the ground state at the vertical link or the horizontal link, we have the
following fusion rules.

$\displaystyle I\to
m_{1,\mathbf{p}-\mathbf{j}_{x}}\otimes\overline{m}_{1,\mathbf{p}}\otimes\overline{e}_{2,\mathbf{r}},\quad
I\to\overline{m}_{1,\mathbf{p}-\mathbf{j}_{y}}\otimes{m}_{1,\mathbf{p}}\otimes
e_{3,\mathbf{r}},$ (63)

which is the concise statement of (60). From these fusion rules, jointly with
(62), we have

${m}_{1,\mathbf{p}}\simeq
m_{1,\mathbf{p}-\mathbf{j}_{x}}\otimes\overline{{\mathfrak{E}}}_{2},\quad{m}_{1,\mathbf{p}}\simeq
m_{1,\mathbf{p}-\mathbf{j}_{y}}\otimes\overline{\mathfrak{E}}_{3}.$ (64)

The iterative use of (64) gives

${m}_{1,\mathbf{p}+K_{x}\mathbf{j}_{x}+K_{y}\mathbf{j}_{y}}\simeq{m}_{1,\mathbf{p}}\otimes\overline{{\mathfrak{E}}}_{2}^{K_{x}}\otimes\overline{\mathfrak{E}}_{3}^{K_{y}}\quad(K_{x},K_{y}\in\mathbb{Z}).$
(65)

In order for (65) to be consistent with the periodic boundary condition, the
multiple of magnetic charges ${m}^{\alpha_{xy}}_{1,\mathbf{p}}$ rather than a
single one are allowed. Here,
$\alpha_{xy}\vcentcolon=\frac{N}{\gcd(N,L_{x},L_{y})}$. One can verify that
the charge ${m}^{\alpha_{xy}}_{1,\mathbf{p}}$ at any coordinate of $\hat{x}$
and $\hat{y}$ is generated by $\overline{{\mathfrak{E}}}_{2}$,
$\overline{\mathfrak{E}}_{3}$, and

${m}^{\alpha_{xy}}_{1,\mathbf{p}_{0}}(\vcentcolon={\mathfrak{M}}_{1})\quad{\rm
with}\ \ \mathbf{p}_{0}=\left(\frac{3}{2},\frac{3}{2}\right),$ (66)

which indicates that the magnetic charge induced by applying the Pauli
operator $X_{1}$ on the ground state is characterized by
$(\mathfrak{M}_{1})^{s_{1}}$ with $s_{1}\in\mathbb{Z}_{\gcd(N,L_{x},L_{y})}$.
Similar discussion leads to that the magnetic charge induced by applying the
Pauli operator $X_{2}$ and $X_{3}$ on the ground state is labeled by

${m}^{\alpha_{x}}_{2,\mathbf{p}_{0}}(\vcentcolon={\mathfrak{M}}_{2})\quad{\rm
and}\quad{m}^{\alpha_{y}}_{3,\mathbf{p}_{0}}(\vcentcolon={\mathfrak{M}}_{3}),$
(67)

respectively, with $\alpha_{x}=\frac{N}{\gcd(N,L_{x})}$ and
$\alpha_{y}=\frac{N}{\gcd(N,L_{y})}$. Overall, there are in total
$\gcd(N,L_{x},L_{y})\times\gcd(N,L_{x})\times\gcd(N,L_{y})$ magnetic charges.

To recap the argument, by investigating the fusion rules of the fractional
excitations, we identify the number of the superselection sectors of the model
on the torus geometry. The excitation is labeled by

$(\mathfrak{E}_{1})^{w_{1}}(\mathfrak{E}_{2})^{w_{2}}(\mathfrak{E}_{3})^{w_{3}}(\mathfrak{M}_{1})^{s_{1}}(\mathfrak{M}_{2})^{s_{2}}(\mathfrak{M}_{3})^{s_{3}}\;\;$
(68)

with $w_{i}\in\mathbb{Z}_{N}\;(i=1,2,3)$,
$s_{1}\in\mathbb{Z}_{\gcd(N,L_{x},L_{y})}$,
$s_{2}\in\mathbb{Z}_{\gcd(N,L_{x})}$ and $s_{3}\in\mathbb{Z}_{\gcd(N,L_{y})}$.
Thus, the number of the superselection sectors, which is equivalent to the GSD
on the torus, is given by

${\rm
GSD}=N^{3}\times\gcd(N,L_{x},L_{y})\times\gcd(N,L_{x})\times\gcd(N,L_{y}).$
(69)

The model that we have considered in this section corresponds to the following
BF theory [32]:

$\mathcal{L}_{{\rm dipole}-x,y}=\frac{N}{2\pi}a\wedge
db+\sum_{I=x,y}\frac{N}{2\pi}A^{I}\wedge dc^{I}+\frac{N}{2\pi}A^{I}\wedge
b\wedge e^{I}.$ (70)

Here, the field $a,b,c^{I},A^{I}$ denotes a $U(1)$ $1$-form gauge field. Also,
$e^{I}$ describes the foliation field, i.e, $e^{x}=dx,e^{y}=dy$. The gauge
field $a$, $A^{x}$, and $A^{y}$ is associated with global and dipole symmetry
in the $x$\- and $y$-direction, respectively. The BF theory (70) shows the
same GSD as (69) [32]. Also, it is interesting to note that a stabilizer model
in a different form obtained by Higgssing the tensor gauge theory [44], which
is a generalized Maxwell theory preserving dipole charges, exhibits the same
GSD as (69) [20, 21].

##  5 Construction of other fracton models

So far we have demonstrated a concrete lattice realization of the topological
phases with multipole symmetries. One can apply the approach discussed in the
previous sections to other fracton model, such as the $\mathbb{Z}_{N}$ exotic
theory and $X$-cube model.

###  5.1 Exotic $\mathbb{Z}_{N}$ theory

To start, let us first introduce BF theory description of the fracton model,
which can be obtained by the BF theory of topological phases with multipole
symmetries with a slight modification. If we replace the gauge field $c^{I}$
with $\phi^{I}e^{I}$ ($I$ is not summed over) in (70), where $\phi^{I}$ is a
$0$-form field, we have

$\mathcal{L}_{EX}=\frac{N}{2\pi}a\wedge
db+\sum_{I=x,y}\frac{N}{2\pi}A^{I}\wedge d\phi^{I}\wedge
e^{I}+\frac{N}{2\pi}A^{I}\wedge b\wedge e^{I},$ (71)

which describes the exotic $\mathbb{Z}_{N}$ gauge theory [40]. Compared with
(70), which has the dipole symmetries, theory (71) possesses an additional
symmetry: the theory is invariant under $A^{I}\to A^{I}+e^{I}g$ with $g$ being
an arbitrary function. This symmetry is known as subsystem symmetry [45, 46,
13], which is of importance in the context of the fracton topological phases.
A theory with such a symmetry admits a symmetry operation on a sub-manifold
rather than entire space. Manifestation of the subsystem symmetry in the
theory (71) is that a gauge invariant operator constructed by the gauge field
$A^{I}$ is mobile only along a sub-manifold, which forbids moving in the
direction parallel to $e^{I}$.

The last term in (71) has the exactly the same form as the one we have
encountered in the previous sections, the foliated DW twist term. Also, the
second term that enters (71) describes foliated topological field theories
each of which corresponds to a spontaneously symmetry breaking phase [47].
Based on this observation, it would be tantalizing to implement the similar
procedure as we did in the previous argument to construct the exotic
$\mathbb{Z}_{N}$ gauge theory via gauging SPT phases.

(a)

(b)

(c)

Figure 4: (a) 2D square lattice and arrays of 1D chains stacked along the
$x$-direction (red squares and dashed lines) and the ones along the
$y$-direction (green triangles and dashed lines). We set the orientation of
the links in the same way as the previous cases. (b) Top: The lattice
configuration after gauging the global $\mathbb{Z}_{N}$ symmetry where gauge
fields reside on links. Bottom: Terms that constitute the Hamiltonian (75).
(c) Examples of the logical operators given in (78).

To this end, we think of a 2D square lattice jointly with arrays of 1D chains
stacked along both of $x$\- and $y$-directions, as demonstrated in Fig. 4(a)
(note the distinction between the configurations in Fig. 4(a) and Fig. 3(a).)
with periodic boundary condition being imposed (i.e.,
$\hat{x}\sim\hat{x}+L_{x}$, $\hat{y}\sim\hat{y}+L_{y}$). We define the Hilbert
space of $\mathbb{Z}_{N}$ qubit on a node of the lattice by $\ket{a}_{i}$ with
$i$ distinguishing the 2D lattice ($i=1$ black dot in Fig. 4(a)) and that of
the 1D chain stacked along $x$\- ($i=2$ red square in Fig. 4(a)) and
$y$-direction ($i=3$ green triangle in Fig. 4(a)). We also introduce
Hamiltonian of this lattice as

$\widetilde{H}_{0}=-\sum_{\mathbf{r}}\left(\tau^{X}_{1,\mathbf{r}}+\tau^{X}_{2,\mathbf{r}}+\tau^{X}_{3,\mathbf{r}}\right).$
(72)

Even though the Hamiltonian has the same form as (53), there is a crucial
difference between the two; in the present case there is no link connecting
between the states $\ket{a}_{2,\mathbf{r}}$ and
$\ket{a}_{2,\mathbf{r}\pm\mathbf{j}_{x}}$ in the horizontal direction.
Likewise, there is no link that connects the states $\ket{a}_{3,\mathbf{r}}$
and $\ket{a}_{3,\mathbf{r}\pm\mathbf{j}_{y}}$ in the vertical direction.

We implement the same CZ gates (54) on the Hamiltonian (72), yielding

$\displaystyle\widetilde{U}^{\dagger}_{CZ}\widetilde{H}_{0}\widetilde{U}_{CZ}=-\sum_{\mathbf{r}}\biggl{[}(\tau^{Z})^{\dagger}_{3,\mathbf{r}-\mathbf{j}_{x}}(\tau^{Z})^{\dagger}_{2,\mathbf{r}-\mathbf{j}_{y}}\tau^{X}_{1,\mathbf{r}}\tau^{Z}_{2,\mathbf{r}}\tau^{Z}_{3,\mathbf{r}}+(\tau^{Z})^{\dagger}_{1,\mathbf{r}+\mathbf{j}_{y}}\tau^{X}_{2,\mathbf{r}}\tau^{Z}_{1,\mathbf{r}}+Z^{\dagger}_{1,\mathbf{r}+\mathbf{j}_{x}}\tau^{X}_{3,\mathbf{r}}\tau^{Z}_{1,\mathbf{r}}\biggr{]}.$
(73)

The model has the one $\mathbb{Z}_{N}$ global and two $\mathbb{Z}_{N}$
subsystem symmetries generated by

$\displaystyle\prod_{\mathbf{r}}\tau^{X}_{1,\mathbf{r}},\quad\prod_{\hat{y}=1}^{L_{y}}\tau^{X}_{2,\mathbf{r}}\;\;(1\leq\hat{x}\leq
L_{x}),\quad\prod_{\hat{x}=1}^{L_{x}}\tau^{X}_{3,\mathbf{r}}\;\;(1\leq\hat{y}\leq
L_{y}).$ (74)

Following the same procedure outlined around (30), one can gauge these
symmetries. The gauged Hamiltonian reads (see also Fig. 4(b))

$\displaystyle H_{EX}=$ $\displaystyle-$
$\displaystyle\biggl{[}\sum_{\mathbf{r}}{V}_{1,\mathbf{r}}+\sum_{\mathbf{p}}{P}_{1,\mathbf{p}}\biggr{]}$
(75) $\displaystyle-$
$\displaystyle\sum_{\hat{x}}H_{fr,x}(\hat{x})-\sum_{\hat{y}}H_{Fr,y}(\hat{y})+h.c.~{}.$

The first line describes the $\mathbb{Z}_{N}$ toric code consisting of the
terms given in (57) and (58) whereas the second line corresponds to
ferromagnetic chains stacked along $x$\- and $y$-direction:

$\displaystyle
H(\hat{x})_{Fr,x}=\sum_{\hat{y}}Z_{1,\mathbf{r}+\mathbf{j}_{y}}X_{2,\mathbf{r}}X^{\dagger}_{2,\mathbf{r}-\mathbf{j}_{y}},\quad
H(\hat{y})_{Fr,y}=\sum_{\hat{x}}Z_{1,\mathbf{r}+\mathbf{j}_{x}}X_{3,\mathbf{r}}X^{\dagger}_{3,\mathbf{r}-\mathbf{j}_{x}}.$
(76)

Note that due to the DW twist terms obtained by gauging the foliated SPT
phases, an additional Pauli operator $Z_{i}$ is attached to the terms (see
also Fig. 4(c)). The ground state of the model (75) is the projected state,
satisfying

$\displaystyle{V}_{1,\mathbf{r}}\ket{\Omega}={P}_{1,\mathbf{p}}\ket{\Omega}=\ket{\Omega}\quad\forall\mathbf{r},\mathbf{p},$
$\displaystyle
Z_{1,\mathbf{r}+\mathbf{j}_{y}}X_{2,\mathbf{r}}X^{\dagger}_{2,\mathbf{r}-\mathbf{j}_{y}}\ket{\Omega}=\ket{\Omega},\quad
Z_{1,\mathbf{r}+\mathbf{j}_{x}}X_{3,\mathbf{r}}X^{\dagger}_{3,\mathbf{r}-\mathbf{j}_{x}}\ket{\Omega}=\ket{\Omega}\quad\forall\mathbf{r}.$
(77)

One can evaluate the GSD of the model on the torus geometry. The model (75)
admits the following logical operators i.e., operators that commute with the
Hamiltonian (75):

$\displaystyle\zeta_{\mathbf{r}}$ $\displaystyle\vcentcolon=$ $\displaystyle
X_{3,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}X^{\dagger}_{2,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}},$
$\displaystyle\eta_{x}(\hat{y})$ $\displaystyle\vcentcolon=$
$\displaystyle\prod_{\hat{x}=1}^{L_{x}}Z_{2,\mathbf{r}-\frac{\mathbf{j}_{y}}{2}}X^{\dagger}_{1,\mathbf{r}-\frac{\mathbf{j}_{y}}{2}}X_{1,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}\quad(1\leq\hat{y}\leq
L_{y}),$ $\displaystyle\;\eta_{y}(\hat{x})$ $\displaystyle\vcentcolon=$
$\displaystyle\prod_{\hat{y}=1}^{L_{y}}Z_{3,\mathbf{r}-\frac{\mathbf{j}_{x}}{2}}X^{\dagger}_{1,\mathbf{r}-\frac{\mathbf{j}_{x}}{2}}X_{1,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}\quad(1\leq\hat{x}\leq
L_{x}).$ (78)

Examples of these operators are shown in Fig. 4(c). The form of the last two
operators resemble the Wilson loop introduced in the exotic $\mathbb{Z}_{N}$
theory constructed by a non-contractible loop of a symmetric tensor gauge
field [40]. There are a few constrains on the operators (78). To see this, we
multiply the term $V_{1,\mathbf{r}}^{\dagger}$ on entire vertex gives

$\displaystyle
1=\prod_{\mathbf{r}}V_{1,\mathbf{r}}^{\dagger}=\prod_{\hat{y}}\eta_{x}(\hat{y})\times\prod_{\hat{x}}\eta_{y}(\hat{x}),$
(79)

where we have used (77). Also, by combination of the terms that enter the
Hamiltonian (75), together with (77) gives

$\displaystyle
1=\zeta_{\mathbf{r}}\zeta^{\dagger}_{\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}\zeta^{\dagger}_{\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}\zeta_{\mathbf{r}+\frac{\mathbf{j}_{x}}{2}+\frac{\mathbf{j}_{y}}{2}}.$
(80)

These constraints (79) and (80) indicate that there are $L_{x}+L_{y}-1$
distinct number of the logical operators $\zeta_{\mathbf{r}}$, and the same
number of distinct logical operators $\eta_{x}(\hat{y})$ and
$\eta_{y}(\hat{x})$. One can show that pairs of such logical operators, each
of which generates $\mathbb{Z}_{N}$ Heisenberg algebra gives the GSD as

${\rm GSD}=N^{L_{x}+L_{y}-1}.$ (81)

This is in agreement with [40].

###  5.2 X-cube model

One can analogously establish the X-cube model in our approach. To this end,
we recall the BF theory description of the X-cube model [34]. Consider the
following BF theory defined in $(3+1)$d

$\mathcal{L}_{{\rm dipole}-xyz}=\frac{N}{2\pi}a\wedge
db+\sum_{I=x,y,z}\frac{N}{2\pi}A^{I}\wedge dc^{I}+\frac{N}{2\pi}A^{I}\wedge
b\wedge e^{I},$ (82)

where the fields $a,A^{I}$ ($b,c^{I}$) represent $U(1)$ $1$-form ($2$-form)
gauge fields and $e^{x}=dx,e^{y}=dy,e^{z}=dz$. This BF theory is the $(3+1)$d
analogue of (70), where the gauge field $a,A^{I}$ is associated with global
and dipole in the $I$th direction, respectively. We replace the $2$-form gauge
field $c^{I}$ with $dB^{I}\wedge e^{I}$ ($I$ is not summed over) where $B^{I}$
denotes a $U(1)$ $1$-form gauge field. With this replacement, the BF theory
(82) becomes

$\mathcal{L}=\frac{N}{2\pi}a\wedge db+\sum_{I=x,y,z}\frac{N}{2\pi}A^{I}\wedge
dB^{I}\wedge e^{I}+\frac{N}{2\pi}A^{I}\wedge b\wedge e^{I},$ (83)

which is known as the effective field theory of the X-cube model [34].
Similarly to (71), the theory has the subsystem symmetry; the theory is
invariant under $A^{I}\to A^{I}+ge^{I}$.

Note that the second term in (83) describes the 2D toric codes, forming
stacked layers and also that the last term in (83) has the same form as the
one we have discussed so far, allowing us to carry out our approach to build
up the UV model which corresponds to the field theory (83). We emphasize that
although realization of the X-cube model via gauging the sub-system symmetry
was discussed previously [13], we believe that it is still useful to present
our way to construct the X-cube model via gauging the foliated SPT phases.
This is because the way we construct the model has a broad perspective in the
sense that it is applicable to both of topological phases with multipole
symmetries and other fracton models, and allows us to make a transparent
understanding of the topological field theory of these phases by interpreting
the coupling term $A^{I}\wedge b\wedge e^{I}$ as the foliated DW twist term.
Furthermore, we believe that our approach may find an application to
realization of the fracton topological phases with more generic group
symmetries including non-Abelian, as our approach complies well with the group
cohomology that describes SPT phases. We leave more details of these issues
for future studies.

(a)

(b)

Figure 5: (a) 3D cubic lattice and layers of 2D lattices. Each node of the 2D
lattice is located at the center of the cube. (b) Side view of the lattice
from the $x$-axis is shown in the middle large panel. The face of the cubic
lattice and the layer of 2D lattice indicated by the yellow arrow in (a)
constitute a 2D lattice. The CZ gates are implemented between states on the
the layer of the 2D lattices and the ones on the face of the cube. Example of
the CZ gates are depicted in the small panels on the left and right. The node
with alphabet corresponds to the one in (a).

With these regards in mind, now we turn to constructing the X-cube model via
gauging foliated SPT phases in the case of $N=2$ with system size $L_{x}\times
L_{y}\times L_{z}$ for simplicity. To do this, we envisage a cubic lattice and
layers of the 2D lattices stacked in the $x$, $y$, and $z$-directions, which
corresponds to blue, red, and green layers in Fig. 5(a) respectively. Each
stacked layers are separated by a unit lattice spacing and each node of the 2D
lattice is located at the center of the cube. We introduce qubits on each link
of the cubic lattice and the ones on each node of the layers of the 2D
lattices. The Hilbert space of the qubit is denoted by $\ket{a}_{i}$ with
$i=0,1,2,3$ and $a=0,1$, where the index $i$ distinguishes the qubits on the
cube ($i=0$) and the ones on the 2D lattices stacked in the $x(i=1)$-,
$y(i=2)$-, and $z(i=3)$-directions. We define the following paramagnet
Hamiltonian:

$\displaystyle
H=-\sum_{\mathbf{r}}\left(\tau^{X}_{0,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}+\tau^{X}_{0,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}+\tau^{X}_{0,\mathbf{r}+\frac{\mathbf{j}_{z}}{2}}\right)-\sum_{\mathbf{r}}\sum_{i=1}^{3}\tau^{X}_{i,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}+\frac{\mathbf{j}_{y}}{2}+\frac{\mathbf{j}_{z}}{2}},$
(84)

where $\mathbf{j}_{z}:=(0,0,1)$,
$\tau^{X}_{0,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}$ is the Pauli operator
acting on the qubit at the link of the cubic lattice with the coordinate
$\mathbf{r}+\frac{\mathbf{j}_{x}}{2}$ and the other Pauli operators are
similarly defined.

We think of implementing the CZ gates (21) on qubits which are located at a
face of the cubic lattice and the ones on the adjacent layer of the 2D
lattice. Instead of writing explicitly the terms, we make use of visual
illustrations given in Figs. 5(a) and 5(b). The CZ gates act on qubits on a 2D
lattice, say the one stacked along the $x$-direction, and the closest face of
the cube in the $+x$-direction, parallel to the $yz$-plane. In Fig. 5(a), they
correspond to the layer and face marked by the yellow arrow. Viewing from the
$x$-axis, the layer and face form a 2D lattice where qubits are placed on each
node and link, as shown in Fig. 5(b). The CZ gates act on a qubit defined on a
node and four qubits that surround it, as well as on a qubit located on a link
and two qubits that are connected with it. The former corresponds to $a\sim e$
in the left of Fig. 5(b) and the latter does to $a$, $d,g$, and $a,e,f$ in the
right of Fig. 5(b). Such manipulation reminds us of the realization of cluster
states of the SPT phase preserving the $0$-form and $1$-form global
symmetries. The CZ gates are implemented on qubits defined on a different
layers stacked along the $y$\- and $z$-directions and the ones on the closest
faces in the similar manner.

The model has the symmetries generated by the following operators

$\displaystyle\prod_{e\perp S}\tau^{X}_{0,e},$
$\displaystyle\prod_{\hat{y},\hat{z}}\tau^{X}_{1,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}+\frac{\mathbf{j}_{y}}{2}+\frac{\mathbf{j}_{z}}{2}}\;(1\leq\hat{x}\leq
L_{x}),\prod_{\hat{z},\hat{x}}\tau^{X}_{2,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}+\frac{\mathbf{j}_{y}}{2}+\frac{\mathbf{j}_{z}}{2}}\;(1\leq\hat{y}\leq
L_{y}),\;\prod_{\hat{x},\hat{y}}\tau^{X}_{3,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}+\frac{\mathbf{j}_{y}}{2}+\frac{\mathbf{j}_{z}}{2}}\;(1\leq\hat{z}\leq
L_{z}).$ (85)

Here the first product represents the multiplication of the $\tau^{X}_{0}$
operators on links that cross a closed surface $S$ in the dual sites, which
corresponds to a global $1$-form symmetry [48, 49] whereas the last three
products do to the subsystem symmetries. One can gauge these symmetries.
Following the same procedure outlined in Sec. 3.1, we gauge the last three
symmetries in (85) on each layer of the 2D lattice.

Regarding gauging the $1$-form $\mathbb{Z}_{2}$ global symmetry, we follow the
analogous steps for gauging ordinary $\mathbb{Z}_{N}$ symmetries. Since the
procedure closely parallels the one in the previous cases [50], we outline how
to carry out gauging succinctly. We introduce $2$-form gauge fields which
reside on each face $f$ of the cube with the Pauli operator $X_{0,f},Z_{0,f}$,
and impose the Gauss law constraint on each link $e$, $G_{e}\ket{\rm
phys}=\ket{\rm phys}$ with

$G_{e}=\tau^{X}_{0,e}\prod_{f|e\subset\partial f}X_{0,f},$

where the product corresponds to multiplication of four Pauli operators that
belong to the link. We further minimally couple the quadratic term on links
that belong to a face, $\prod_{e\subset\partial f}\tau^{Z}_{0,e}$ to the
$2$-form gauge field via $Z_{0,f}\prod_{e\subset\partial f}\tau^{Z}_{0,e}$ and
impose the flux-less condition of the 2-form gauge fields.

(a)

(b)

(c)

Figure 6: (a) Configurations of the qubits after gauging and redefining the
lattice grid. (b) The terms defined in (87) and (88) that constitute the
Hamiltonian (LABEL:lattice_fl). (c) Configurations of loops given in (92).

After gauging the symmetries (85), we also redefine the lattice grid to make
the figure more visually-friendly. We do so in such a way that gauge fields
are located on each link of the newly-defined cubic lattice; the gauge fields
introduced by implementing gauging the first symmetry in (85) is located on
links of the cubic lattice (black dots in Fig. 6(a)) whereas the ones
introduced by the second, third, and fourth symmetry in (85) distinguished by
the index $i=1$, $i=2$ and $i=3$, is placed on links belonging to the face of
the cubic lattice, in the $yz$-, $zx$-, and $xy$-planes, respectively (blue,
red, and green dots in Fig. 6(a)).

We finally arrive at the following gauged Hamiltonian:

$\displaystyle H_{{\rm gauged},Xc}=$ $\displaystyle-$
$\displaystyle\biggl{[}\sum_{\mathbf{r}}V_{0,\mathbf{r}}+\sum_{\hat{z}}\sum_{\mathbf{p}_{xy}(\hat{z})}B_{0,\mathbf{p}_{xy}(\hat{z})}+\sum_{\hat{x}}\sum_{\mathbf{p}_{yz}(\hat{x})}B_{0,\mathbf{p}_{yz}(\hat{x})}+\sum_{\hat{y}}\sum_{\mathbf{p}_{zx}(\hat{y})}B_{0,\mathbf{p}_{zx}(\hat{y})}\biggr{]}$
$\displaystyle-$
$\displaystyle\sum_{\hat{x}}\biggl{[}\sum_{\hat{y},\hat{z}}V_{1,\mathbf{r}}+\sum_{\mathbf{p}_{yz}(\hat{x})}B_{1,\mathbf{p}_{yz}(\hat{x})}\biggr{]}-\sum_{\hat{y}}\biggl{[}\sum_{\hat{x},\hat{z}}V_{2,\mathbf{r}}+\sum_{\mathbf{p}_{zx}(\hat{y})}B_{2,\mathbf{p}_{zx}(\hat{y})}\biggr{]}-\sum_{\hat{z}}\biggl{[}\sum_{\hat{x},\hat{y}}V_{3,\mathbf{r}}+\sum_{\mathbf{p}_{xy}(\hat{z})}B_{3,\mathbf{p}_{zx}(\hat{z})}\biggr{]}.$

The individual term of the first line in (LABEL:lattice_fl) is defined by

$\displaystyle V_{0,\mathbf{r}}$ $\displaystyle\vcentcolon=$
$\displaystyle\prod_{t=\pm
1}Z_{0,\mathbf{r}+t\frac{\mathbf{j}_{x}}{2}}Z_{0,\mathbf{r}+t\frac{\mathbf{j}_{y}}{2}}Z_{0,\mathbf{r}+t\frac{\mathbf{j}_{z}}{2}},$
$\displaystyle B_{0,\mathbf{p}_{xy}(\hat{z})}$ $\displaystyle\vcentcolon=$
$\displaystyle
Z_{1,\mathbf{p}_{xy}(\hat{z})-\frac{\mathbf{j}_{x}}{2}}Z_{2,\mathbf{p}_{xy}(\hat{z})-\frac{\mathbf{j}_{y}}{2}}\prod_{t=\pm
1}X_{0,\mathbf{p}_{xy}(\hat{z})+t\frac{\mathbf{j}_{x}}{2}}X_{0,\mathbf{p}_{xy}(\hat{z})+t\frac{\mathbf{j}_{y}}{2}},$
$\displaystyle B_{0,\mathbf{p}_{yz}(\hat{x})}$ $\displaystyle\vcentcolon=$
$\displaystyle
Z_{2,\mathbf{p}_{yz}(\hat{x})-\frac{\mathbf{j}_{z}}{2}}Z_{3,\mathbf{p}_{yz}(\hat{x})-\frac{\mathbf{j}_{y}}{2}}\prod_{t=\pm
1}X_{0,\mathbf{p}_{yz}(\hat{x})+t\frac{\mathbf{j}_{y}}{2}}X_{0,\mathbf{p}_{yz}(\hat{x})+t\frac{\mathbf{j}_{z}}{2}},$
$\displaystyle B_{0,\mathbf{p}_{zx}(\hat{y})}$ $\displaystyle\vcentcolon=$
$\displaystyle
Z_{3,\mathbf{p}_{zx}(\hat{y})-\frac{\mathbf{j}_{x}}{2}}Z_{1,\mathbf{p}_{zx}(\hat{y})-\frac{\mathbf{j}_{z}}{2}}\prod_{t=\pm
1}X_{0,\mathbf{p}_{zx}(\hat{y})+t\frac{\mathbf{j}_{x}}{2}}X_{0,\mathbf{p}_{zx}(\hat{y})+t\frac{\mathbf{j}_{z}}{2}},$
(87)

which are depicted in the first four terms in Fig. 6(b). Here, we have defined
the coordinate of a plaquette in the $xy$, $yz$ and $zx$-plane by
$\mathbf{p}_{xy}(\hat{z})=(\hat{x}+\frac{1}{2},\hat{y}+\frac{1}{2},\hat{z})$,
$\mathbf{p}_{yz}(\hat{x})=(\hat{x},\hat{y}+\frac{1}{2},\hat{z}+\frac{1}{2})$
and
$\mathbf{p}_{zx}(\hat{y})=(\hat{x}+\frac{1}{2},\hat{y},\hat{z}+\frac{1}{2})$.
Also, the individual term of the second, third and fourth line in
(LABEL:lattice_fl) reads

$\displaystyle
V_{1,\mathbf{r}}\vcentcolon=Z_{0,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}\prod_{t=\pm
1}X_{1,\mathbf{r}+t\frac{\mathbf{j}_{y}}{2}}X_{1,\mathbf{r}+t\frac{\mathbf{j}_{z}}{2}},$
$\displaystyle\quad B_{1,\mathbf{p}_{yz}(\hat{x})}\vcentcolon=\prod_{s=\pm
1}Z_{1,\mathbf{p}_{yz}(\hat{x})+s\frac{\mathbf{j}_{y}}{2}}Z_{1,\mathbf{p}_{yz}(\hat{x})+s\frac{\mathbf{j}_{z}}{2}},$
$\displaystyle
V_{2,\mathbf{r}}\vcentcolon=Z_{0,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}\prod_{t=\pm
1}X_{2,\mathbf{r}+t\frac{\mathbf{j}_{x}}{2}}X_{2,\mathbf{r}+t\frac{\mathbf{j}_{z}}{2}},$
$\displaystyle\quad B_{2,\mathbf{p}_{zx}(\hat{y})}\vcentcolon=\prod_{s=\pm
1}Z_{2,\mathbf{p}_{zx}(\hat{y})+s\frac{\mathbf{j}_{x}}{2}}Z_{2,\mathbf{p}_{zx}(\hat{y})+s\frac{\mathbf{j}_{z}}{2}},$
$\displaystyle
V_{3,\mathbf{r}}\vcentcolon=Z_{0,\mathbf{r}\frac{\mathbf{j}_{z}}{2}}\prod_{t=\pm
1}X_{3,\mathbf{r}+t\frac{\mathbf{j}_{x}}{2}}X_{3,\mathbf{r}+t\frac{\mathbf{j}_{z}}{2}},$
$\displaystyle\quad B_{3,\mathbf{p}_{xy}(\hat{z})}\vcentcolon=\prod_{s=\pm
1}Z_{3,\mathbf{p}_{xy}(\hat{z})+s\frac{\mathbf{j}_{x}}{2}}Z_{3,\mathbf{p}_{xy}(\hat{z})+s\frac{\mathbf{j}_{y}}{2}}.$
(88)

These terms are portrayed in the last four terms in Fig. 6(b).

The ground state of the Hamiltonian (LABEL:lattice_fl) is the projected state,
satisfying that all of the eigenvalues of the terms in (87) and (88) are
trivial, e.g., $V_{0,\mathbf{r}}=1$. Essentially, the first line in
(LABEL:lattice_fl) describes the 3D toric code whereas the terms in the second
line in (LABEL:lattice_fl) do the layers of the 2D toric codes stacked along
the $x$-, $y$-, and $z$-directions. The crucial difference between the toric
codes and the model (LABEL:lattice_fl) is that a few Pauli operators are
multiplied with the plaquette and vertex terms of the toric codes, which alter
the statistics of the fractional excitations, mirroring the fact that we
accommodate the foliated DW twist terms.

We investigate the GSD of the model (LABEL:lattice_fl) on a 3D torus with
system size $L_{x}\times L_{y}\times L_{z}$ by counting the number of the
distinct non-contractible loops of the electric charges. Practically, we
evaluate the number of the distinct operators of the Pauli $Z_{i}$ operators
that form the non-contractible loops. The model admits a number of non-
contractible loops constructed by the product of the Pauli $Z_{i}$ operators:

$\displaystyle
W_{1e,y}(\hat{x})=\prod_{\hat{y}=1}^{L_{y}}Z_{1,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}},$
$\displaystyle\quad
W_{1e,z}(\hat{x})=\prod_{\hat{z}=1}^{L_{z}}Z_{1,\mathbf{r}+\frac{\mathbf{j}_{z}}{2}}\qquad(1\leq\hat{x}\leq
L_{x}),$ $\displaystyle
W_{2e,z}(\hat{y})=\prod_{\hat{z}=1}^{L_{z}}Z_{2,\mathbf{r}+\frac{\mathbf{j}_{z}}{2}},$
$\displaystyle\quad
W_{2e,x}(\hat{y})=\prod_{\hat{x}=1}^{L_{x}}Z_{2,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}}\qquad(1\leq\hat{y}\leq
L_{y}),$ $\displaystyle
W_{3e,x}(\hat{z})=\prod_{\hat{x}=1}^{L_{x}}Z_{3,\mathbf{r}+\frac{\mathbf{j}_{x}}{2}},$
$\displaystyle\quad
W_{3e,y}(\hat{z})=\prod_{\hat{y}=1}^{L_{y}}Z_{1,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}\qquad(1\leq\hat{z}\leq
L_{z}).$ (89)

Examples are shown in the left of Fig. 6(c). Naively there are
${2L_{x}+2L_{y}+2L_{z}}$ distinct loops, however, there are a few constraints
on them, reducing the number.777 One might wonder that the model admits non-
contractible membrane operators which are found in the 3D toric code, such as
$W_{0e,xy}=\prod_{\hat{x}=1}^{L_{x}}\prod_{\hat{y}=1}^{L_{y}}Z_{0,\mathbf{r}+\frac{\mathbf{j}_{z}}{2}}$.
However, such membranes become trivial as multiplication of the terms in (88)
gives identity. Indeed, we have
$1=\prod_{\hat{y}=1}^{L_{y}}\prod_{\hat{z}=1}^{L_{z}}V_{1,\mathbf{r}}=W_{0e,xy}$.
Indeed, multiplication of the operators in (87) gives

$\displaystyle 1$ $\displaystyle=$
$\displaystyle\prod_{\hat{x},\hat{y}}B_{0,\mathbf{p}_{xy}(\hat{z})}=\prod_{\hat{y}}W_{2e,x}(\hat{y})\times\prod_{\hat{x}}W_{1e,y}(\hat{x}),$
$\displaystyle 1$ $\displaystyle=$
$\displaystyle\prod_{\hat{y},\hat{z}}B_{0,\mathbf{p}_{yz}(\hat{x})}=\prod_{\hat{z}}W_{3e,y}(\hat{z})\times\prod_{\hat{y}}W_{2e,z}(\hat{y}),$
$\displaystyle 1$ $\displaystyle=$
$\displaystyle\prod_{\hat{z},\hat{x}}B_{0,\mathbf{p}_{zx}(\hat{y})}=\prod_{\hat{z}}W_{3e,x}(\hat{z})\times\prod_{\hat{x}}W_{1e,z}(\hat{x}),$
(90)

where we have used the fact that the ground states satisfies that the
eigenvalues of the operators (87) are trivial. Such constraints come from the
fact that a few Pauli operators are attached to the terms in the 2D toric
codes, due to the foliated DW twist terms that we introduced. Hence, the
distinct number of non-contractible loops of the $\mathbb{Z}_{2}$ electric
charges is given by ${2L_{x}+2L_{y}+2L_{z}}-3$, leading to

${\rm GSD}=2^{2L_{x}+2L_{y}+2L_{z}-3},$ (91)

which agrees with the GSD found in [13].

One can analogously count the number of distinct loops of the magnetic
charges, i.e., the number of non-contractible loops constructed by $X_{i}$’s.
Examples of such loops are

$\displaystyle
W_{2m,x}(\hat{y})=\prod_{\hat{x}}X_{2,\mathbf{r}+\mathbf{j}_{y}+\frac{\mathbf{j}_{z}}{2}}X_{2,\mathbf{r}+\frac{\mathbf{j}_{z}}{2}}Z_{0,\mathbf{r}+\mathbf{j}_{y}+\frac{\mathbf{j}_{z}}{2}},\quad
W_{2m,3m,x}(\hat{y})=\prod_{\hat{x}}X_{2,\mathbf{r}+\frac{\mathbf{j}_{z}}{2}}X_{3,\mathbf{r}+\frac{\mathbf{j}_{y}}{2}}$
(92)

which are depicted in Fig. 6(c). We can find other loops in the similar form.
By the similar argument around (90), it can be verified that there are
$2L_{x}+2L_{y}+2L_{z}-3$ distinct loops. By using these loops as well as the
ones we have considered (89), one can make $2L_{x}+2L_{y}+2L_{z}-3$ pairs of
loops each of which generates the $\mathbb{Z}_{2}$ Heisenberg algebra,
arriving at the GSD (91).

##  6 Discussion and conclusion

Originally proposed in the context of quantum information [11, 12], the
fracton topological phases have now diverse research interests, involving
various branches of physics, such as condensed matter physics, and high-energy
physics. Motivated by the ultimate goal to construct a complete theoretical
framework to describe such phases, the concept of symmetries has been
developed. Focusing on one of new symmetries, multipole symmetry, we have
demonstrated a way to build up topological phases which respect such a
symmetry. Take-home message in this work is that via gauging SPT phases with
global $\mathbb{Z}_{N}$ symmetries, they become the foliated DW twist terms
which yield unconventional fractional excitations found in the multipole
topological phases. Furthermore, this construction gives a clear physical
understanding of the term $A^{I}\wedge b\wedge e^{I}$ in a field theory
description of the multipole topological orders. The way we construct phases
can also be implemented in other fracton models. Before closing this section,
we make a few remarks and future research directions, regarding the present
work.

###  6.1 Comments on entanglement entropy

Recently, several works studied entanglement entropy of topological phases
with multipole symmetries. It is widely known that the entanglement entropy of
a contractible disk geometry $A$ in a conventional $(2+1)$d topologically
ordered phase reads $S_{A}=A_{\rm area}-\gamma$ [51, 52]. While the first term
is the so-called area law term, proportional to the perimeter of the disk, the
second sub-leading order term is particularly important. It is known as the
topological entanglement entropy which relates to a topological order via
$\gamma=\log\sqrt{\sum_{a}d_{a}^{2}}$, where $d_{a}$ denotes a quantum
dimension of an anyon labeled by $a$. It was shown that the topological
entanglement entropy of a contractible disk geometry in a $\mathbb{Z}_{N}$
topological phase with multipole symmetries gives $\gamma=m\log N$
($m\in\mathbb{Z}$), which is in contrast with the fact that the number of
superselection sectors, that is, the number of distinct fractional excitations
depends on the greatest common divisor between $N$ and the system size [53,
54].

Such a result can be understood by the form of the BF theory of topological
phases with multipole symmetries (7) and (70) together with the interpretation
of the coupling term that we make throughout this paper, i.e., interpretation
of the term $A^{I}\wedge b\wedge e^{I}$ as the foliated DW twist terms.
Generally, it was shown that the DW twist terms alter the fractional
statistics of anyons but do not change the total quantum dimension. Indeed, it
was proven that in a concrete lattice model of a quantum double with DW twist
terms the topological entanglement entropy of a disk geometry takes the same
values as the one in the case without the DW twist [55]. By this result,
together with our observation that the BF theory of the topological phases
with multipole symmetries consists of layers of the toric codes with the
foliated DW twist terms, leads to that the topological entanglement entropy of
such phases has the form $\gamma=k\log N$, where $k$ is the number of layers
of the toric codes (more explicitly, the number of degrees of freedom of the
multipoles such as global charge, dipole, quadrupole, e.t.c.).

As an example, the topological entanglement entropy of a disk in the UV
stabilizer models respecting the dipole symmetry in the $x$-direction was
found to be $\gamma=2\log N$ [53], which is consistent with the form of the BF
theory (7) as the theory consists of the two copies of the toric codes with
the DW twist; regardless of the presence of the DW twist, the total quantum
dimension of the excitations amounts to be the one of the two copies of the
toric codes, hence $\gamma=2\log N$.

###  6.2 Future directions

There are several future directions regarding the present work. Studying
topological phases with multipole symmetries in view of quantum information,
such as quantum error correction [8], is a widely open problem. Especially, it
could be theoretically and practically important to address how stable the
model considered in this work is against errors e.g., bit or phase flips and
how the multipole symmetry plays a role compared with the conventional
topological stabilizer codes. Investigating observable in the context of
quantum information other than entanglement entropy, such as quantum
negativity [56], could be an another interesting direction.

It would be interesting to extend our analysis to the case of arrays of SPT
phases with general group symmetry $G$ to realize more general foliated DW
twist terms, as our approach nicely fits into the group cohomology.
Especially, in the case of $G$ being non-Abelian, the investigation on the
compatibility condition with the periodic boundary condition becomes
complicated, which could yield richer multipole or fracton topological phases.

It would also be intriguing to incorporate other types of symmetries into our
construction. The gauge transformation of the multipole symmetries (4)
resembles the one found in $2$-group symmetry [57] in that different gauge
fields are related with one another via gauge transformation. Studying whether
one can construct topological phases with higher group symmetry, or even
richer topological phases where both of higher-group and multipole symmetries
are incorporated by our approach would contribute to making better
understanding of new topological phases with various symmetries.

Recently it was found that several field theories related to fractons have ’t
Hooft anomalies and their anomaly inflows [58, 59, 60, 61, 62, 38]. It would
be illuminating to study relations among them and foliated BF theories. Last
but not least, some foliated BF theories themselves may also have ’t Hooft
anomalies. It would be also interesting to study this aspect in future.

## Acknowledgement

We thank Qiang Jia, Ken Shiozaki and Shimamori Souichiro for helpful
discussion. H. E. is supported by KAKENHI-PROJECT-23H01097. M. H. is supported
by MEXT Q-LEAP, JSPS Grant-in-Aid for Transformative Research Areas (A)
“Extreme Universe" JP21H05190 [D01] and JSPS KAKENHI Grant Number 22H01222. M.
H. and T. N. are supported by JST PRESTO Grant Number JPMJPR2117. T. N is
supported by JST, the establishment of university fellowships towards the
creation of science technology innovation, Grant Number JPMJFS2123.

## References

  * [1] D. C. Tsui, H. L. Stormer, and A. C. Gossard, “Two-dimensional magnetotransport in the extreme quantum limit,” Phys. Rev. Lett., vol. 48, pp. 1559–1562, 1982.
  * [2] V. Kalmeyer and R. B. Laughlin, “Equivalence of the resonating-valence-bond and fractional quantum hall states,” Phys. Rev. Lett., vol. 59, pp. 2095–2098, Nov 1987.
  * [3] X. G. Wen, F. Wilczek, and A. Zee, “Chiral spin states and superconductivity,” Phys. Rev. B, vol. 39, pp. 11413–11423, Jun 1989.
  * [4] N. Read and S. Sachdev, “Large-n expansion for frustrated quantum antiferromagnets,” Physical review letters, vol. 66, no. 13, p. 1773, 1991.
  * [5] J. M. Leinaas and J. Myrheim, “On the theory of identical particles,” Il Nuovo Cimento B (1971-1996), vol. 37, no. 1, pp. 1–23, 1977.
  * [6] R. B. Laughlin, “Anomalous quantum hall effect: An incompressible quantum fluid with fractionally charged excitations,” Phys. Rev. Lett., vol. 50, pp. 1395–1398, May 1983.
  * [7] A. Kitaev, “Fault-tolerant quantum computation by anyons,” Annals of Physics, vol. 303, no. 1, pp. 2–30, 2003.
  * [8] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, “Topological quantum memory,” Journal of Mathematical Physics, vol. 43, no. 9, pp. 4452–4505, 2002.
  * [9] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma Rev. Mod. Phys., vol. 80, pp. 1083–1159, 2008.
  * [10] R. S. K. Mong et al., “Universal Topological Quantum Computation from a Superconductor-Abelian Quantum Hall Heterostructure,” Phys. Rev. X, vol. 4, no. 1, p. 011036, 2014.
  * [11] C. Chamon, “Quantum glassiness in strongly correlated clean systems: An example of topological overprotection,” Phys. Rev. Lett., vol. 94, p. 040402, Jan 2005.
  * [12] J. Haah, “Local stabilizer codes in three dimensions without string logical operators,” Phys. Rev. A, vol. 83, p. 042330, Apr 2011.
  * [13] S. Vijay, J. Haah, and L. Fu, “Fracton topological order, generalized lattice gauge theory, and duality,” Phys. Rev. B, vol. 94, p. 235157, Dec 2016.
  * [14] T. Griffin, K. T. Grosvenor, P. Hořava, and Z. Yan, “Scalar field theories with polynomial shift symmetries,” Communications in Mathematical Physics, vol. 340, pp. 985–1048, 2015.
  * [15] M. Pretko, “The Fracton Gauge Principle,” Phys. Rev. B, vol. 98, no. 11, p. 115134, 2018.
  * [16] A. Gromov, “Towards classification of fracton phases: The multipole algebra,” Phys. Rev. X, vol. 9, p. 031035, Aug 2019.
  * [17] H. Ma, M. Hermele, and X. Chen, “Fracton topological order from the higgs and partial-confinement mechanisms of rank-two gauge theory,” Phys. Rev. B, vol. 98, p. 035111, Jul 2018.
  * [18] P. Gorantla, H. T. Lam, N. Seiberg, and S.-H. Shao, “Global dipole symmetry, compact lifshitz theory, tensor gauge theory, and fractons,” Phys. Rev. B, vol. 106, p. 045112, Jul 2022.
  * [19] A. Jain and K. Jensen, “Fractons in curved space,” SciPost Phys., vol. 12, no. 4, p. 142, 2022.
  * [20] S. D. Pace and X.-G. Wen, “Position-dependent excitations and uv/ir mixing in the ${\mathbb{z}}_{N}$ rank-2 toric code and its low-energy effective field theory,” Phys. Rev. B, vol. 106, p. 045145, 2022.
  * [21] Y.-T. Oh, J. Kim, E.-G. Moon, and J. H. Han, “Rank-2 toric code in two dimensions,” Physical Review B, vol. 105, no. 4, p. 045128, 2022.
  * [22] H. Ebisu and B. Han, “Anisotropic higher rank $\mathbb{Z}_{N}$ topological phases on graphs,” SciPost Phys., vol. 14, p. 106, 2023.
  * [23] H. Ebisu, “Symmetric higher rank topological phases on generic graphs,” Phys. Rev. B, vol. 107, p. 125154, Mar 2023.
  * [24] X. Chen, Z.-C. Gu, Z.-X. Liu, and X.-G. Wen, “Symmetry-protected topological orders in interacting bosonic systems,” Science, vol. 338, no. 6114, pp. 1604–1606, 2012.
  * [25] X. Chen, Z.-C. Gu, Z.-X. Liu, and X.-G. Wen, “Symmetry protected topological orders and the group cohomology of their symmetry group,” Phys. Rev. B, vol. 87, p. 155114, Apr 2013.
  * [26] X. Chen, Y.-M. Lu, and A. Vishwanath, “Symmetry-protected topological phases from decorated domain walls,” Nature communications, vol. 5, no. 1, p. 3507, 2014.
  * [27] R. Dijkgraaf and E. Witten, “Topological gauge theories and group cohomology,” Communications in Mathematical Physics, vol. 129, no. 2, pp. 393 – 429, 1990.
  * [28] M. d. W. Propitius, “Topological interactions in broken gauge theories,” arXiv preprint hep-th/9511195, 1995.
  * [29] M. Levin and Z.-C. Gu, “Braiding statistics approach to symmetry-protected topological phases,” Phys. Rev. B, vol. 86, p. 115109, Sep 2012.
  * [30] J. C. Wang, Z.-C. Gu, and X.-G. Wen, “Field-theory representation of gauge-gravity symmetry-protected topological invariants, group cohomology, and beyond,” Phys. Rev. Lett., vol. 114, p. 031601, Jan 2015.
  * [31] C. Wang and M. Levin, “Topological invariants for gauge theories and symmetry-protected topological phases,” Phys. Rev. B, vol. 91, p. 165119, Apr 2015.
  * [32] H. Ebisu, M. Honda, and T. Nakanishi, “Foliated BF theories and Multipole symmetries,” arXiv e-prints, p. arXiv:2310.06701, Oct. 2023.
  * [33] D. Aasen, D. Bulmash, A. Prem, K. Slagle, and D. J. Williamson, “Topological defect networks for fractons of all types,” Phys. Rev. Res., vol. 2, p. 043165, Oct 2020.
  * [34] K. Slagle, D. Aasen, and D. Williamson, “Foliated field theory and string-membrane-net condensation picture of fracton order,” SciPost Phys., vol. 6, p. 043, 2019.
  * [35] K. Slagle, “Foliated quantum field theory of fracton order,” Phys. Rev. Lett., vol. 126, p. 101603, Mar 2021.
  * [36] K. Ohmori and S. Shimamura, “Foliated-exotic duality in fractonic BF theories,” SciPost Phys., vol. 14, no. 6, p. 164, 2023.
  * [37] R. C. Spieler, “Exotic field theories for (hybrid) fracton phases from imposing constraints in foliated field theory,” JHEP, vol. 09, p. 178, 2023.
  * [38] W. Cao and Q. Jia, “Symmetry TFT for Subsystem Symmetry,” arXiv preprint arXiv:2310.01474, 2023.
  * [39] L. Bhardwaj, S. Schäfer-Nameki, and A. Tiwari, “Unifying constructions of non-invertible symmetries,” SciPost Phys., vol. 15, p. 122, 2023.
  * [40] N. Seiberg and S.-H. Shao, “Exotic symmetries, duality, and fractons in 2+ 1-dimensional quantum field theory,” SciPost Physics, vol. 10, no. 2, p. 027, 2021.
  * [41] Y. Hirono, M. You, S. Angus, and G. Y. Cho, “A symmetry principle for gauge theories with fractons,” arXiv preprint arXiv:2207.00854, 2022.
  * [42] P. Gorantla, H. T. Lam, N. Seiberg, and S.-H. Shao, “A modified Villain formulation of fractons and other exotic theories,” J. Math. Phys., vol. 62, no. 10, p. 102301, 2021.
  * [43] B. Yoshida, “Gapped boundaries, group cohomology and fault-tolerant logical gates,” Annals of Physics, vol. 377, pp. 387–413, 2017.
  * [44] M. Pretko, “Generalized electromagnetism of subdimensional particles: A spin liquid story,” Phys. Rev. B, vol. 96, p. 035119, Jul 2017.
  * [45] A. Paramekanti, L. Balents, and M. P. A. Fisher, “Ring exchange, the exciton bose liquid, and bosonization in two dimensions,” Phys. Rev. B, vol. 66, p. 054526, Aug 2002.
  * [46] N. Seiberg, “Field Theories With a Vector Global Symmetry,” SciPost Phys., vol. 8, no. 4, p. 050, 2020.
  * [47] A. Kapustin and N. Seiberg, “Coupling a qft to a tqft and duality,” Journal of High Energy Physics, vol. 2014, no. 4, pp. 1–45, 2014.
  * [48] R. Raussendorf, S. Bravyi, and J. Harrington, “Long-range quantum entanglement in noisy cluster states,” Physical Review A, vol. 71, no. 6, p. 062313, 2005.
  * [49] D. Gaiotto, A. Kapustin, N. Seiberg, and B. Willett, “Generalized global symmetries,” Journal of High Energy Physics, vol. 2015, no. 2, pp. 1–62, 2015.
  * [50] A. Kapustin and R. Thorngren, “Higher symmetry and gapped phases of gauge theories,” Algebra, Geometry, and Physics in the 21st Century: Kontsevich Festschrift, pp. 177–202, 2017.
  * [51] M. Levin and X.-G. Wen, “Detecting topological order in a ground state wave function,” Phys. Rev. Lett., vol. 96, p. 110405, Mar 2006.
  * [52] A. Kitaev and J. Preskill, “Topological entanglement entropy,” Phys. Rev. Lett., vol. 96, p. 110404, Mar 2006.
  * [53] H. Ebisu, “Entanglement entropy of higher rank topological phases,” arXiv:2302.11468 (2023).
  * [54] J. Kim, Y.-T. Oh, D. Bulmash, and J. H. Han, “Unveiling UV/IR Mixing via Symmetry Defects: A View from Topological Entanglement Entropy,” arXiv e-prints, p. arXiv:2310.09425, Oct. 2023.
  * [55] X. Wen, H. He, A. Tiwari, Y. Zheng, and P. Ye, “Entanglement entropy for (3+1)-dimensional topological order with excitations,” Phys. Rev. B, vol. 97, p. 085147, Feb 2018.
  * [56] Y. A. Lee and G. Vidal, “Entanglement negativity and topological order,” Phys. Rev. A, vol. 88, p. 042318, Oct 2013.
  * [57] C. Córdova, T. T. Dumitrescu, and K. Intriligator, “Exploring 2-Group Global Symmetries,” JHEP, vol. 02, p. 184, 2019.
  * [58] F. J. Burnell, T. Devakul, P. Gorantla, H. T. Lam, and S.-H. Shao, “Anomaly inflow for subsystem symmetries,” Phys. Rev. B, vol. 106, no. 8, p. 085113, 2022.
  * [59] S. Yamaguchi, “Gapless edge modes in (4+1)-dimensional topologically massive tensor gauge theory and anomaly inflow for subsystem symmetry,” PTEP, vol. 2022, no. 3, p. 033B08, 2022.
  * [60] M. Honda and T. Nakanishi, “Scalar, fermionic and supersymmetric field theories with subsystem symmetries in d + 1 dimensions,” JHEP, vol. 03, p. 188, 2023.
  * [61] W. Cao, L. Li, M. Yamazaki, and Y. Zheng, “Subsystem non-invertible symmetry operators and defects,” SciPost Phys., vol. 15, no. 4, p. 155, 2023.
  * [62] X. Huang, “A Chern-Simons theory for dipole symmetry,” SciPost Phys., vol. 15, no. 4, p. 153, 2023.