# Subdimensional criticality: condensation of lineons and planons in the
X-cube model

Ethan Lake  Department of Physics, Massachusetts Institute of Technology,
Cambridge, MA, 02139, USA  Michael Hermele  Department of Physics and Center
for Theory of Quantum Matter, University of Colorado, Boulder, CO, 80309, USA

###### Abstract

We study quantum phase transitions out of the fracton ordered phase of the
$\mathbb{Z}_{N}$ X-cube model. These phase transitions occur when various
types of sub-dimensional excitations and their composites are condensed. The
condensed phases are either trivial paramagnets, or are built from stacks of
$d=2$ or $d=3$ deconfined gauge theories, where $d$ is the spatial dimension.
The nature of the phase transitions depends on the excitations being
condensed. Upon condensing dipolar bound states of fractons or lineons, for
$N\geq 4$ we find stable critical points described by decoupled stacks of
$d=2$ conformal field theories. Upon condensing lineon excitations, when $N>4$
we find a gapless phase intermediate between the X-cube and condensed phases,
described as an array of $d=1$ conformal field theories. In all these cases,
effective subsystem symmetries arise from the mobility constraints on the
excitations of the X-cube phase and play an important role in the analysis of
the phase transitions.

##  I Introduction and summary

Fracton phasesChamon (2005); Haah (2011); Yoshida (2013); Vijay et al. (2016);
Pretko (2017); Nandkishore and Hermele (2019); Pretko et al. (2020a) are an
exciting class of states of matter in three dimensions which defy many
expectations for the types of behavior one expects to see in zero-temperature
quantum systems. They have subextensively large ground state degeneracy,
excitations whose motions are restricted to lie along low-dimensional
submanifolds of space, and provide examples of gapped phases not described
within the context of conventional topological quantum field theory.111We will
not attempt a review of the many interesting recent advances in the field of
fractons; for an introduction to the literature the reader may consult the
reviewsNandkishore and Hermele (2019); Pretko et al. (2020a) and the
references therein.

Fracton phases can be separated into two types, depending on the nature of
their excitations.Vijay et al. (2016) Type-I fracton phases possess mobile
excitations, which are typically _lineons_ or _planons_ that can move along
one- or two-dimensional submanifolds of space, respectively. One paradigmatic
example of a type-I fracton phase is the X-Cube (XC) model.Vijay et al. (2016)
Type-II phases by contrast do not have any non-trivial mobile excitations,
with Haah’s codeHaah (2011) being the classic example.

An important question in fracton physics is where fracton phases may arise in
the phase diagrams of various systems. One way of approaching this problem is
to examine what sorts of continuous quantum phase transitions can occur
between a phase with fracton order and a more conventional phase of matter.
This question has already been partially addressed in the literature. Previous
works have identified first-order quantum phase transitions out of fracton
ordered phases via a range of techniques,Devakul et al. (2018); Slagle and Kim
(2017a); Mühlhauser et al. (2020) and studied the breakdown of fracton order
in series expansions.Poon and Liu (2020) Ref. Vijay, 2017 exploited a duality
to 3d $\mathbb{Z}_{N}$ gauge theory to propose a continuous transition between
the $\mathbb{Z}_{N}$ XC phase and a trivial gapped phase; we comment on this
proposal in Sec. VI. (Here and throughout the paper, we refer to spatial
dimension rather than space-time dimension, unless explicitly stated
otherwise.) A recent work introduced the notion of hybrid fracton orders,
which combine the fully mobile excitations of 3d topological orders with the
restricted-mobility excitations of fracton phases, and found continuous phase
transitions between such hybrid phases and fracton-ordered
phases.Tantivasadakarn et al. (2021) However, fracton order is present on both
sides of these phase transitions. Other recent works have studied continuous
quantum phase transitions in two-dimensional systems with subsystem symmetry
but without fracton order.You et al. (2020a, 2021); Zhou et al. (2021)

In this work, our focus will be on continuous phase transitions out of fracton
orders, where the continuous nature and stability of the phase transitions may
be demonstrated rigorously. We will study both critical points, which can be
accessed by tuning a single parameter, as well as stable gapless phases
proximate to fracton orders. Given that the excitations in a fracton phase
have restricted mobility, it is natural to imagine that any putative critical
point would necessarily involve a rather unconventional field theory
description, perhaps of a similar flavor to the novel types of field theories
that have appeared when analyzing fracton phases themselves.Slagle and Kim
(2017b); Slagle et al. (2019); You et al. (2020b); Seiberg and Shao (2020a, b,
c); Gorantla et al. (2020); Slagle (2021); Hsin and Slagle (2021)

The easiest place to start when thinking about continuous phase transitions
out of fracton phases is to examine phase transitions in type-I fracton phases
where certain types of excitations with restricted mobility are condensed. In
these types of phase transitions, it is however not clear whether or not one
should expect a critical point characterized by scale invariance, and to what
extent the standard theory of critical phenomena should be applicable.

In this paper we will make some first steps towards addressing these
questions. We will focus on what appear to be the simplest examples of
continuous phase transitions out of fracton orders, viz. various types of
condensation transitions in the $\mathbb{Z}_{N}$ X-Cube (XC) model.Vijay et
al. (2016); Vijay (2017) The nature of the condensed phase varies depending on
which type of excitation is condensed. In some of the phase transitions, the
condensed phase is trivial, while in others it possesses topological or
fracton order (distinct from that of the $\mathbb{Z}_{N}$ X-cube phase).

Figure 1:  A schematic of the phase transitions we consider. Tuning a
parameter $g$ drives a quantum phase transition between the $\mathbb{Z}_{N}$
X-Cube model and a condensate of non-trivial lineon or planon excitations. The
resulting condensed phase is either trivial or is some form of deconfined
gauge theory. A critical point or intermediate gapless phase can be described
as a stack of 2d or 1d CFTs.

Taking inspiration from the many approaches which construct fracton phases by
coupling various well-understood lower-dimensional phases together,Ma et al.
(2017); Vijay (2017); Halász et al. (2017); Prem et al. (2019); Aasen et al.
(2020); Wen (2020); Sullivan et al. (2020); Slagle et al. (2019) our approach
to studying phase transitions in these models will be to look for critical
points, or intermediate gapless phases, which can be understood within the
context of 1d or 2d critical phenomena. These critical points and gapless
phases arise when condensing planon or lineon excitations of the
$\mathbb{Z}_{N}$ XC model, and can be described as stacks / arrays of familiar
2d / 1d conformal field theories (see figure 1 for a schematic).

An interesting feature of the critical points and gapless phases we study is
the emergence of effective subsystem symmetries. By definition, subsystem
symmetries are those which act on degrees of freedom lying only in a subspace,
which may be a line, a plane, or even a fractal. Subsystem symmetries are
interesting in their own right and also play a variety of important roles in
the theory of fracton phases. The effective subsystem symmetries in the
theories we study arise from the gauge structure of the $\mathbb{Z}_{N}$
X-cube phase, and are robust to arbitrary perturbations, even though such
symmetries are not assumed to be present microscopically.

sector | condensate | condensed phase
---|---|---
lineon | $\mathfrak{e}^{z}$ |  $\mathbb{Z}_{N,2}$ stack
| $\mathfrak{e}^{y,z}$ | trivial
| $\mathfrak{d}^{z}_{\mathfrak{e}}$ | two $\mathbb{Z}_{N,2}$ stacks
| $\mathfrak{d}^{y,z}_{\mathfrak{e}}$ |  $\mathbb{Z}_{N,3}$ \+ $\mathbb{Z}_{N,2}$ stack
| $\mathfrak{d}^{x,y,z}_{\mathfrak{e}}$ | $\mathbb{Z}_{N,3}^{2}$
fracton | $\mathfrak{d}^{z}_{\mathfrak{m}}$ | anisotropic model
| $\mathfrak{d}^{y,z}_{\mathfrak{m}}$ |  $\mathbb{Z}_{N,2}$ stack
| $\mathfrak{d}^{x,y,z}_{\mathfrak{m}}$ | $\mathbb{Z}_{N,3}$
Table 1:  The phases obtained when condensing various types of excitations in
the $\mathbb{Z}_{N}$ XC model. $\mathfrak{e}^{a}$ denotes a lineon mobile in
the $a$-direction, while $\mathfrak{d}^{a}_{\mathfrak{e}}$ and
$\mathfrak{d}^{a}_{\mathfrak{m}}$ are dipolar bound states of lineons and
fractons, respectively, with dipole moment in the $a$-direction; see Sec. II
for more details. $\mathbb{Z}_{N,d}$ denotes a deconfined $\mathbb{Z}_{N}$
gauge theory in $d$ spatial dimensions, a “stack” refers to a decoupled set of
theories, and “trivial” refers to a paramagnet lacking any topological order.
We have identified continuous phase transitions involving lineon condensation
provided that $N>4$, and continous transitions involving dipole condensation
provided that $N>3$. “Anisotropic model” refers to the model of lineons and
planeons discussed in Ref. Shirley et al., 2019.

A summary of this paper is as follows. We start in section II by reviewing
some facts about the $\mathbb{Z}_{N}$ XC model which will be relevant for the
following sections. In section III we study what happens when dipoles of
fractons, which are planons, are condensed. In order to identify the resulting
condensed phases we employ a dual representation of the XC model in terms of a
generalized gauge theory, which is detailed in section III.1. Condensing
fracton dipoles allows single fractons, which in the XC model are immobile, to
move freely in directions parallel to the moments of the condensed dipoles.
When dipoles with a single direction of dipole moment are condensed the
resulting phase is an anisotropic model studied in Ref. Shirley et al., 2019,
when two orthogonal directions of dipoles are condensed one finds a stack of
deconfined 2d $\mathbb{Z}_{N}$ gauge theories, and when all dipoles are
condensed the resulting phase is a deconfined $\mathbb{Z}_{N}$ gauge theory,
with single fractons serving as the gauge charges. In section III.3 we examine
the nature of the condensation transition. When $N<4$ the condensation
transitions are likely to be generically first order, while when $N\geq 4$
continuous phase transitions exist, and are described by stacks of critical 2d
XY models.

Section IV is analogous to the previous section, except we instead consider
the condensation of lineon dipoles. When lineon dipoles condense the mobility
restrictions on single lineons are relaxed, and single lineons become capable
of moving freely in two or three dimensions. When a single orientation of
lineon dipole condenses the resulting phase is two interpenetrating stacks of
2d $\mathbb{Z}_{N}$ gauge theories, when two orientations of dipoles condense
one obtains a 3d $\mathbb{Z}_{N}$ gauge theory and a stack of 2d
$\mathbb{Z}_{N}$ gauge theories, and when all dipoles condense the result is a
3d $\mathbb{Z}_{N}^{2}$ gauge theory. These results are summarized in table 1,
and are obtained by utilizing a generalized gauge theory which makes the
nature of the condensed phases very explicit, and which is explained in detail
in section IV.1. The condensation phase transitions in this case are exactly
the same as those that occur when condensing fracton dipoles.

In section V we discuss the slightly simpler problem of condensing lineons.
When only one type of lineon condenses the resulting phase is a stack of 2d
$\mathbb{Z}_{N}$ gauge theories, while if more condense one obtains a trivial
paramagnet. In section V.1 we introduce yet another generalized gauge theory,
which is used to obtain a field theory description of lineon condensation. In
sections V.2 and V.3 we show that for $N\leq 4$ the condensation transition is
likely to always be first order, while for $N>4$ it is possible for an
intermediate gapless phase to exist in between the X-Cube and condensed
phases. This gapless phase is described by an array of strongly coupled
Luttinger liquids, and is related to the 3d sliding phases discussed in Ref.
Mukhopadhyay et al., 2001. The stability of this gapless phase is somewhat
delicate, and we argue that (1) by adding inter-chain derivative couplings the
phase can be made stable down to very low temperatures and (2) a truly stable
phase with more complicated inter-chain couplings should exist in principle.

Finally, section VI contains a conclusion and a brief discussion of some
remaining open questions.

Figure 2:  A graphical representation of the terms in the Hamiltonian of the
$\mathbb{Z}_{N}$ XC model.

##  II Preliminaries

Before starting in earnest, we will quickly summarize a few important
properties of the $\mathbb{Z}_{N}$ X-cube (XC) model,Vijay et al. (2016);
Vijay (2017) as well as some of the notation that we will be making frequent
use of. $Z,X$ will denote $\mathbb{Z}_{N}$ clock and shift operators, which
obey the relations

$ZX=e^{{2\pi}i/N}XZ,\qquad Z^{\dagger}X=e^{-{2\pi}i/N}XZ^{\dagger},$ (1)

where $N$ will be inferred from context. When discussing microscopic
Hamiltonians, we will be making use of various dual representations. These
dual representations will be in terms of discrete $\mathbb{Z}_{N}$ matter
spins and discrete $\mathbb{Z}_{N}$ gauge fields. The clock and shift
operators for the matter fields will always be written in serif font
($\mathsf{z},\mathsf{x}$), while those for the gauge fields will always be in
mathcal ($\mathcal{Z},\mathcal{X}$).

The Hilbert space of the $\mathbb{Z}_{N}$ XC model is built from
$\mathbb{Z}_{N}$ spins placed on the links of a 3d cubic lattice. We will
write the Hamiltonian for the commuting projector limit of the XC model as

$H_{XC}=-g\sum_{i}(A_{i}^{x}+A_{i}^{y}+A^{z}_{i})-K\sum_{c}B_{c}+h.c.$ (2)

where the sums over $i$ and $c$ run over the sites and minimal cubes of the
lattice, respectively. The operators $A^{a}_{i}$ are defined as products of
$Z$ operators over the four links touching a given vertex $i$, while the
$B_{c}$ operators are defined as products of $X$ operators over the twelve
links of the cube $c$. Half of the $X$ and $Z$ operators appearing in the
$A^{a}_{i}$ and $B_{c}$ terms are Hermitian conjugated in the manner shown in
figure 2, which ensures that all of the terms in (2) commute, and that

$A^{x}_{i}A^{y}_{i}A^{z}_{i}=\mathbf{1}$ (3)

at every site $i$.

Sometimes it is convenient to refer to links of the lattice in an oriented
manner. We imagine that each link comes with a given orientation in the $+x$,
$+y$ or $+z$ direction. Then, denoting by $\langle ij\rangle$ a pair of
nearest-neighbor sites, and letting ${\mathbf{\hat{a}}}$ be the unit vector
pointing from $i$ to $j$, we define $X_{\langle ij\rangle}=X$ if
${\mathbf{\hat{a}}}$ is parallel to the link’s given orientation, and
$X_{\langle ij\rangle}=X^{\dagger}$ if ${\mathbf{\hat{a}}}$ is anti-parallel
to the given orientation. Note that $X_{\langle
ij\rangle}=X^{\dagger}_{\langle ji\rangle}$. The same notation and conventions
are used to define $Z_{{\langle ij\rangle}}$. For simplicity, we will use $l$
to denote links when keeping track of the orientation is not necessary.

The excitations of the Hamiltonian (2) are divided into two sectors. The first
sector consists of lineon excitations and their composites. Single lineons
correspond to violations of two of the $A^{a}_{i}$ terms, and come in three
species, which we will denote as
$\mathfrak{e}^{x},\mathfrak{e}^{y},\mathfrak{e}^{z}$. An $\mathfrak{e}^{x}$
excitation at a given site $i$ corresponds to a violation of the two terms
$A^{y}_{i}$ and $A^{z}_{i}$, and similarly for
$\mathfrak{e}^{y},\mathfrak{e}^{z}$. Our convention is that acting on a ground
state with $X_{\langle ij\rangle}$, where _e.g._ $\langle ij\rangle$ is
directed along the positive or negative $x$-direction, creates an
$\mathfrak{e}_{i}^{x}$ lineon at site $i$ and a
$(\mathfrak{e}_{j}^{x})^{\dagger}$ antilineon at site $j$. Each isolated
$\mathfrak{e}^{a}$ lineon may only move along the ${\mathbf{\hat{a}}}$
direction; any attempts to move it along the directions normal to
${\mathbf{\hat{a}}}$ necessarily result in the creation of additional
excitations. String operators for the lineons are formed from products of
$X_{\langle ij\rangle}$ operators in their direction of motion. Due to the
relation (3), the lineons obey the fusion rule

$\mathfrak{e}_{i}^{x}\mathfrak{e}^{y}_{i}\mathfrak{e}^{z}_{i}=1.$ (4)

Because of this fusion rule, a lineon $\mathfrak{e}^{x}$ moving along the
${\mathbf{\hat{x}}}$ direction may split into an
$(\mathfrak{e}^{y})^{\dagger}$ antilineon and an
$(\mathfrak{e}^{z})^{\dagger}$ antilineon, and likewise for the other two
species.

The mobility of individual and multiple lineons can be understood by viewing
$\mathfrak{e}^{a}$ as carrying a unit $\mathbb{Z}_{N}$ charge in each of the
two lattice planes that intersect at the $a$-axis along which
$\mathfrak{e}^{a}$ is free to move. (This was discussed in Ref. Pai and
Hermele, 2019 for the $\mathbb{Z}_{2}$ X-cube model; the generalization to
$\mathbb{Z}_{N}$ is straightforward.) The $\mathbb{Z}_{N}$ charge in each
plane is conserved in the sense that it cannot be changed by any local
operator. This gives rise to effective planar $\mathbb{Z}_{N}$ subsystem
symmetries that will be important in understanding critical points where
lineons and their composites condense.

Another interesting feature of the XC model is the behavior of dipolar222Since
lineons are $\mathbb{Z}_{N}$ objects, the “dipole moment” of two lineons is
rather ill-defined. In what follows we will however abuse terminology and
refer to a bound state of two lineons on neighboring vertices as a lineon
“dipole”. bound states of lineons. Consider a bound state of a pair of
$\mathfrak{e}^{b}$ and $(\mathfrak{e}^{b})^{\dagger}$ lineons separated by a
given distance in the ${\mathbf{\hat{a}}}$ direction, where $a\neq b$. Because
the two lineons cannot move in the ${\mathbf{\hat{a}}}$ direction, they cannot
move together and annihilate to vacuum, and indeed it can be shown that an
isolated dipole of this sort is stable, in the sense that it cannot be
destroyed by any local process.333By contrast, if lineon and antilineon are
separated in the ${\mathbf{\hat{b}}}$ direction, they can annihilate, and such
an excitation is trivial. One consequence of the fusion rule (4) is that such
dipolar bound states may move freely within the plane perpendicular to
${\mathbf{\hat{a}}}$, as they may change the direction of their motion within
the plane by having their constituent charges exchange an $\mathfrak{e}^{a}$
lineon. However motion parallel to ${\mathbf{\hat{a}}}$ is not possible, and
consequently such dipolar bound states are planon excitations.

Figure 3:  An example of a string operator which moves a
$\mathfrak{d}^{z}_{\mathfrak{e}}$ lineon dipole (the lattice is not shown for
ease of visualization). The pattern of Hermitian conjugates can be determined
by comparing to the pattern appearing in the definition of the $B_{c}$
operator.

We will denote a “minimal-strength” lineon dipole by
$\mathfrak{d}^{a}_{\mathfrak{e}}$, where ${\mathbf{\hat{a}}}$ is the normal
vector to the plane of motion. Each $\mathfrak{d}^{a}_{\mathfrak{e}}$ is built
of an $\mathfrak{e}^{b}$ and an $(\mathfrak{e}^{b})^{\dagger}$ separated by
one lattice spacing in the ${\mathbf{\hat{a}}}$ direction, where $b\neq a$.
There are no local processes which turn a $\mathfrak{d}^{a}_{\mathfrak{e}}$
dipole into a $\mathfrak{d}^{b\neq a}_{\mathfrak{e}}$ dipole. The string
operators which move lineon dipoles in their planes of motion are “wireframe
operators”, and are given by products of $X_{\langle ij\rangle}$ operators
along the edges of rectuangular prism-like shapes. The best way to understand
the geometry of these string operators is with a picture – see figure 3 for an
example. In terms of planar $\mathbb{Z}_{N}$ charges as described above,
$\mathfrak{d}^{a}_{\mathfrak{e}}$ carries $\mathbb{Z}_{N}$ charges $+1$ and
$-1$ in the two planes normal to ${\mathbf{\hat{a}}}$ in which the constituent
lineon excitations reside.

Figure 4:  A membrane operator supported on the green membrane $\mathcal{M}$,
which is formed as a product of $Z_{\langle ij\rangle}$ operators on the links
marked by the purple circles. Acting on a ground state, this operator creates
four fractons living on the shaded cubes. Blue cubes have charge $+1$, and red
cubes have charge $-1$ (mod $N$).

The second sector of excitations consists of fractons and their composites.
Isolated fractons arise as violations of a single $B_{c}$ term in $H_{XC}$,
and are created with $Z_{\langle ij\rangle}$ operators. Acting with
$Z_{\langle ij\rangle}$ on a ground state creates a configuration of four
fractons at the cubes that contain ${\langle ij\rangle}$ as an edge, with the
$\mathbb{Z}_{N}$ charges of the fractons arranged in a “quadrupolar” pattern.
More generally, acting with a product of $Z_{\langle ij\rangle}$ operators
over all links ${\langle ij\rangle}$ dual to a rectangular membrane
$\mathcal{M}$ of dual-lattice plaquettes creates one fracton at each of the
four corners of the membrane $\mathcal{M}$, as shown in figure 4. A membrane
operator supported on $\mathcal{M}$ obeys $\mathbb{Z}_{N}$ commutation
relations with lineon string operators that intersect $\mathcal{M}$; this fact
enables certain types of statistical braiding-like processes to be defined
between fractons and lineons.Pai and Hermele (2019)

A single isolated fracton cannot move without creating additional excitations.
In more detail, and more generally, the mobility of any collection of fractons
can be understood by thinking of each fracton as carrying a conserved
$\mathbb{Z}_{N}$ charge in the three planes normal to ${\mathbf{\hat{x}}}$,
${\mathbf{\hat{y}}}$ and ${\mathbf{\hat{z}}}$ that intersect at its
position.Vijay et al. (2016); Pai and Hermele (2019) As for lineons, local
operators carry vanishing planar $\mathbb{Z}_{N}$ charges. It is important to
note that the planar charges carried by lineons and fractons are distinct;
these are two independent sectors of excitations.

Figure 5:  An example of a string operator which moves a
$\mathfrak{d}^{z}_{\mathfrak{e}}$ fracton dipole from one location to another.

Similarly to the lineon sector, dipolar444As with lineons, we will continue to
abuse terminology by referring to these bound states as “dipoles”. bound
states of fractons can freely move in the plane normal to their dipole moment.
More specifically, by a fracton dipole we mean a pair of fracton and anti-
fracton excitations separated along the $x$-, $y$\- or $z$-axis. The
separation between the two constituent excitations cannot change by a local
process, and moreover the dipole cannot move along the direction of its
moment. We will denote a “minimal-strength” fracton dipole capable of moving
in a plane normal to ${\mathbf{\hat{a}}}$ as
$\mathfrak{d}^{a}_{\mathfrak{m}}$. The constituent fracton and anti-fracton
excitations in $\mathfrak{d}^{a}_{\mathfrak{m}}$ are separated by a single
lattice spacing along the ${\mathbf{\hat{a}}}$ direction. There is no local
process which can convert a $\mathfrak{d}^{a}_{\mathfrak{m}}$ dipole to a
$\mathfrak{d}^{b\neq a}_{\mathfrak{m}}$ dipole. Fracton dipoles are moved
using products of $Z_{\langle ij\rangle}$ operators along links lying within
the plane of their motion, with an example shown in figure 5. Moreover, just
as for lineon dipoles, $\mathfrak{d}^{a}_{\mathfrak{m}}$ carries $+1$ and $-1$
planar $\mathbb{Z}_{N}$ charges in the planes normal to ${\mathbf{\hat{a}}}$
that contain the constituent fracton and anti-fracton excitations. As a
consequence, the mobility restrictions on fracton and lineon dipoles – if no
other types of excitations are present – are _exactly_ identical.

Due to the relation between the XC model and coupled layers of
$\mathbb{Z}_{N}$ gauge theories,Ma et al. (2017); Vijay (2017) we will refer
to lineon excitations and their composites as “electric” excitations, and
fracton excitations and their composites as “magnetic” excitations.

##  III Fracton dipole condensation

We now turn to discussing the condensation of fracton dipoles in the
$\mathbb{Z}_{N}$ XC model. We will first identify the phases that result upon
condensing the fracton dipoles, using a rather explicit generalized gauge
theory construction with some similarities to the membrane-net model of Ref.
Slagle et al., 2019. Later in section III.3 we will examine the critical
points that occur at the condensation transition.

###  III.1 Condensed phases and dual gauge theories

Before getting into details, let us ask on general grounds what types of
phases we expect to obtain upon condensing various types of fracton dipoles,
while leaving single fractons uncondensed.

Figure 6:  An example of a string operator which moves an isolated fracton
along the ${\mathbf{\hat{x}}}$ direction in the presence of a condensate of
$\mathfrak{d}^{x}_{\mathfrak{m}}$ dipoles.

Consider first condensing all three types of fracton dipoles. Since fracton
dipoles have nontrivial statistical phases with lineons which move in their
plane of motion, all of the lineons will be confined in the condensed phase.
Nevertheless, the resulting phase is not trivial. Indeed, single isolated
fractons remain as well-defined deconfined excitations, as there remains no
local process which can create an isolated fracton. Furthermore by absorbing
dipoles from the condensate, isolated fractons may now move freely in all
three dimensions, via processes like the one shown in figure 6. It is
therefore natural to identify the resulting phase with deconfined 3d
$\mathbb{Z}_{N}$ gauge theory.555In fact, in this case standard 3d
$\mathbb{Z}_{N}$ gauge theory is the only possibility—the set of 3d
topological orders with bosonic gauge charges which have fusion rules given by
a group $G$ are fully enumerated by 3d Dijkgraaf-Witten $G$-gauge theories
(aka twisted $G$-gauge theories),Lan et al. (2018) which in the case of
$G=\mathbb{Z}_{N}$ are classified by
$H^{4}(\mathbb{Z}_{N};U(1))=H^{5}(\mathbb{Z}_{N};\mathbb{Z})=\mathbb{Z}_{1}$,
meaning that $\mathbb{Z}_{N}$ gauge theory is the only possibility.

We can also contemplate condensing only a certain subset of the
$\mathfrak{d}_{\mathfrak{m}}^{a}$ dipoles. First, consider condensing just
$\mathfrak{d}_{\mathfrak{m}}^{z}$. The condensate allows isolated fractons to
move in the ${\mathbf{\hat{z}}}$ direction, so that they become lineons. The
condensate confines $\mathfrak{e}^{x}$ and $\mathfrak{e}^{y}$ lineons, while
$\mathfrak{e}^{z}$ lineons remain deconfined. As before, the
$\mathfrak{e}^{z}$ lineons can pair into $\mathfrak{d}^{x,y}_{\mathfrak{e}}$
dipoles, which are planons.

Next, consider condensing both $\mathfrak{d}^{x}_{\mathfrak{m}}$ and
$\mathfrak{d}^{y}_{\mathfrak{m}}$. This condensate leaves
$\mathfrak{d}^{z}_{\mathfrak{e}}$ dipoles deconfined, but confines the rest of
the lineon sector. In the fracton sector, isolated fractons can move in both
the ${\mathbf{\hat{x}}}$ and ${\mathbf{\hat{y}}}$ directions, but remain
immobile in the ${\mathbf{\hat{z}}}$ direction. The deconfined excitations are
thus single fractons and $\mathfrak{d}^{z}_{\mathfrak{e}}$ dipoles, both of
which are capable of moving only within planes normal to ${\mathbf{\hat{z}}}$.
Since both of these excitations are bosons, the excitation content of the
phase obtained this way is the same as a stack of 2d deconfined
$\mathbb{Z}_{N}$ gauge theories, with the fractons corresponding to the $e$
excitations and the $\mathfrak{d}^{z}_{\mathfrak{e}}$ dipoles corresponding to
the $m$ excitations. We will see later that this is indeed the correct
identification of the condensed phase.

Figure 7:  Conventions for indexing the links on plaquettes. All plaquettes
with a given unit normal are labeled identically.

Now we will show how these expectations can be demonstrated explicitly. For
simplicity of notation, in what follows we will specialize to the $N=2$ case.
The case of general $N$ is similar, and is presented in appendix A. We will
first focus on the case where all species of fracton dipoles condense.

Since in the condensed phase fractons become deconfined mobile particles, it
is helpful to reformulate the XC model as a generalized gauge theory, where
fractons are the gauge charges. To do this it is convenient to perform a
duality mapping on the XC model, with the dual representation being formulated
in terms of spins living on the dual lattice. In what follows, we will use
sans serif symbols to denote elements of the dual lattice, with
$\mathsf{i},\mathsf{l},\mathsf{p},$ and $\mathsf{c}$ standing for dual
vertices, links, plaquettes, and cubes, respectively. Our duality mapping
involves matter qubits $\mathsf{x}_{\mathsf{i}},\mathsf{z}_{\mathsf{i}}$
living on the vertices of the dual lattice, and gauge qubits
$\mathcal{X}_{\mathsf{p}},\mathcal{Z}_{\mathsf{p}}$ living on the plaquettes
of the dual lattice. Sometimes it will be convenient to use $\mathsf{l}$ to
refer to a unit vector along the given link, and $\mathsf{p}$ to refer to a
unit vector normal to the given plaquette, so that for instance
$\mathsf{l}{\,\parallel\,}{\mathbf{\hat{z}}}$ denotes a link parallel to
${\mathbf{\hat{z}}}$, while $\mathsf{p}{\,\parallel\,}{\mathbf{\hat{z}}}$
denotes a plaquette with normal in the ${\mathbf{\hat{z}}}$ direction.

The duality works in much the same way as that between the $\mathbb{Z}_{2}$
toric code and $\mathbb{Z}_{2}$ gauge theory coupled to Ising matter,Kitaev
(2003) with the operators mapping as

$X_{l}\rightarrow\mathcal{X}_{\mathsf{p}},\qquad
Z_{l}\rightarrow\mathcal{Z}_{\mathsf{p}}\prod_{\mathsf{i}\in\mathsf{p}}\mathsf{z}_{\mathsf{i}},$
(5)

where $\mathsf{p}$ is the dual lattice plaquette dual to $l$, and the product
is over the four corners of $\mathsf{p}$. In the dual formulation, there is a
Gauss’ law constraint at each vertex $\mathsf{i}$

$\mathsf{x}_{\mathsf{i}}=\prod_{\mathsf{p}\ni\mathsf{i}}\mathcal{X}_{\mathsf{p}},$
(6)

where the product is over the twelve dual plaquettes meeting $\mathsf{i}$ at a
corner. The Hamiltonian in this representation is

$H_{XC}^{\prime}=-K\sum_{\mathsf{i}}\mathsf{x}_{\mathsf{i}}-g\sum_{\mathsf{c},a}A^{a}_{\mathsf{c}}.$
(7)

Here $A^{a}_{\mathsf{c}}$ is a product of $\mathcal{Z}_{\mathsf{p}}$ over four
of the six plaquettes in the cube $\mathsf{c}$, excluding the two plaquettes
normal to the ${\mathbf{\hat{a}}}$ direction.

We will now modify this model to include a separate field for fracton dipoles.
As such dipoles naturally live on the links of the dual lattice, we introduce
a qubit on each link with Pauli operators $X_{\mathsf{l}},Z_{\mathsf{l}}$ (not
to be confused with the Pauli operators in the original formulation of the XC
model, which live on the direct lattice), where $X_{\mathsf{l}}$ is the dipole
number operator and $Z_{\mathsf{l}}$ creates a fracton dipole at the link
$\mathsf{l}$.

Since single fractons carry gauge charge, a dipole at $\mathsf{l}$ should
carry gauge charge on the vertices at the two ends of $\mathsf{l}$. We
therefore modify the Gauss law constraint to be

$\mathsf{x}_{\mathsf{i}}\prod_{\mathsf{l}\ni
i}X_{\mathsf{l}}=\prod_{\mathsf{p}\ni\mathsf{i}}\mathcal{X}_{\mathsf{p}},$ (8)

where the products are over all the links and plaquettes neighboring the
vertex $\mathsf{i}$. With this Gauss law two fractons at adjacent vertices can
be converted into a dipole via the gauge-invariant operator
$\mathsf{z}_{\mathsf{i}}Z_{\langle\mathsf{i}\mathsf{j}\rangle}\mathsf{z}_{\mathsf{j}}$,
and two dipoles with parallel dipole moments can be created on opposite edges
$\mathsf{l},\mathsf{l}^{\prime}$ of a plaquette $\mathsf{p}$ via the gauge-
invariant operator
$Z_{\mathsf{l}}\mathcal{Z}_{\mathsf{p}}Z_{\mathsf{l}^{\prime}}$.

To analyze what happens when we condense the dipoles, it will be notationally
helpful to fix a standard set of links $\mathsf{l}_{1},\dots,\mathsf{l}_{4}$
for each plaquette of a given orientation. Our conventions for each plaquette
are shown in 7, and are such that for any plaquette $\mathsf{l}_{1}$ is always
parallel to $\mathsf{l}_{2}$, and $\mathsf{l}_{3}$ is always parallel to
$\mathsf{l}_{4}$.

We now consider the following Hamiltonian:

$\displaystyle H_{con}$
$\displaystyle=-K\sum_{\mathsf{i}}\mathsf{x}_{\mathsf{i}}-g\sum_{\mathsf{i},a}A^{a}_{\mathsf{i}}$
(9)
$\displaystyle-h\sum_{\mathsf{l}}X_{\mathsf{l}}-\lambda\sum_{\mathsf{p}}(Z_{\mathsf{l}_{1}}\mathcal{Z}_{\mathsf{p}}Z_{\mathsf{l}_{2}}+Z_{\mathsf{l}_{3}}\mathcal{Z}_{\mathsf{p}}Z_{\mathsf{l}_{4}}).$

The last term proportional to $\lambda$ is a kinetic energy term for the
dipoles, and allows dipoles on links parallel to ${\mathbf{\hat{a}}}$ to move
in planes normal to ${\mathbf{\hat{a}}}$.

When $h/\lambda\gg 1$, we are in the regime where both the dipoles and
fractons are gapped, and the gauge field $\mathcal{Z}$ is deconfined. In this
limit we are clearly in the same phase as $H_{XC}$. In the following we will
instead be interested in the limit $h/\lambda\ll 1$, as in this limit the
dipoles are condensed. In this limit the second term in (9) can be dropped
without changing the low-energy physics. This is because $A^{a}_{\mathsf{i}}$
can be obtained by taking products of the dipole kinetic terms, so that if
these terms are set to $1$ then so are the $A^{a}_{\mathsf{i}}$ terms.

We will now show that in the phase where dipoles are condensed, the
Hamiltonian (9) can be transformed into the standard Hamiltonian of
$\mathbb{Z}_{2}$ gauge theory, plus an unimportant decoupled paramagnet
arising from the plaquette gauge fields. For simplicity and to work at a
solvable point within the dipole-condensed phase, we set $h=g=0$. Note that if
we ignore the plaquette variables $\mathcal{X}_{\mathsf{p}}$, the Gauss law
(8) becomes the ordinary $\mathbb{Z}_{2}$ gauge theory Gauss law constraint.
Our strategy will therefore be to arrive at the $\mathbb{Z}_{2}$ gauge theory
by performing a unitary transformation designed to eliminate the
$\mathcal{X}_{\mathsf{p}}$ variables appearing in (8).

To begin, for each plaquette, define the CNOT-like operators

$U_{\mathsf{p}}\equiv\left(\frac{1+Z_{\mathsf{l}_{1}}Z_{\mathsf{l}_{2}}}{2}\right)+\left(\frac{1-Z_{\mathsf{l}_{1}}Z_{\mathsf{l}_{2}}}{2}\right)X_{\mathsf{p}}.$
(10)

Note that $U_{\mathsf{p}}=U^{\dagger}_{\mathsf{p}}$ is unitary, and is not
gauge invariant — thus we can use it to transform both the Hamiltonian and the
gauge constraint.

The unitary $U_{\mathsf{p}}$ conjugates the various Pauli operators as

$\displaystyle U_{\mathsf{p}}\mathcal{Z}_{\mathsf{p}}U_{\mathsf{p}}$
$\displaystyle=\mathcal{Z}_{\mathsf{p}}Z_{\mathsf{l}_{1}}Z_{\mathsf{l}_{2}}$
(11) $\displaystyle U_{\mathsf{p}}X_{\mathsf{l}_{1,2}}U_{\mathsf{p}}$
$\displaystyle=X_{\mathsf{l}_{1,2}}\mathcal{X}_{\mathsf{p}}$ $\displaystyle
U_{\mathsf{p}}X_{\mathsf{l}_{3,4}}U_{\mathsf{p}}$
$\displaystyle=X_{\mathsf{l}_{3,4}}\text{,}$

with other conjugations trivial.

Noting that $[U_{\mathsf{p}},U_{\mathsf{p}^{\prime}}]=0$, we now perform a
unitary conjugation by the operator

$\mathcal{U}\equiv\prod_{\mathsf{p}}U_{\mathsf{p}}.$ (12)

While the support of $U_{\mathsf{p}}$ operators on adjacent plaquettes does
overlap, $\mathcal{U}$ can be written as a depth-two quantum circuit by
dividing the plaquettes into two subsets. The transformed Hamiltonian then
takes the form

$\displaystyle H^{\prime}_{con}$ $\displaystyle=\mathcal{U}H_{con}\mathcal{U}$
(13)
$\displaystyle=-K\sum_{\mathsf{i}}\mathsf{x}_{\mathsf{i}}-\lambda\sum_{\mathsf{p}}\mathcal{Z}_{\mathsf{p}}-\lambda\sum_{\mathsf{p}}\mathcal{Z}_{\mathsf{p}}\prod_{\mathsf{l}\in\mathsf{p}}Z_{\mathsf{l}}$

Here the first term proportional to $\lambda$ comes from conjugating terms
which hop dipoles between $\mathsf{l}_{1}$ and $\mathsf{l}_{2}$ links, while
the second comes from terms which hop dipoles between $\mathsf{l}_{3}$ and
$\mathsf{l}_{4}$ links. The product in the final term is over the four links
in the perimeter of $\mathsf{p}$.

To compute the transformation of the gauge constraint, note that for a given
vertex $\mathsf{i}$, for each of the twelve plaquettes meeting $\mathsf{i}$,
one of the links $\mathsf{l}_{1},\mathsf{l}_{2}$ has $\mathsf{i}$ as an
endpoint. When the $X_{\mathsf{l}}$ operator on this link is conjugated by
$U$, it produces by (11) an $\mathcal{X}_{\mathsf{p}}$ operator, which cancels
one of the $\mathcal{X}_{\mathsf{p}}$’s appearing in the Gauss law. The
transformed gauge constraint is therefore

$\mathsf{x}_{\mathsf{i}}=\prod_{l\ni i}X_{\mathsf{l}},$ (14)

which is obtained by conjugating both sides of (8) by $\mathcal{U}$. This is
exactly the gauge constraint of $\mathbb{Z}_{2}$ gauge theory.

It only remains to deform $H^{\prime}_{con}$ to the standard solvable
Hamiltonian for the deconfined phase of $\mathbb{Z}_{2}$ gauge theory. This
can obviously be done without closing the gap simply by freezing out the now
gauge-invariant $\mathcal{Z}_{\mathsf{p}}$ degrees of freedom in (13), which
are gapped and have no dynamics. More precisely, we may introduce a parameter
$s\in[0,1]$ and write

$\displaystyle H^{\prime}_{con}(s)$
$\displaystyle\equiv-K\sum_{\mathsf{i}}\mathsf{x}_{\mathsf{i}}-\lambda\sum_{\mathsf{p}}\mathcal{Z}_{\mathsf{p}}$
(15)
$\displaystyle\qquad-\lambda\sum_{\mathsf{p}}[(1-s)\mathcal{Z}_{\mathsf{p}}+s]\prod_{\mathsf{l}\in\mathsf{p}}Z_{\mathsf{l}}.$

Clearly at $s=0$ we have $H^{\prime}_{con}(0)=H^{\prime}_{con}$, the ground
state does not change as a function of $s$, and the gap remains open. At $s=1$
we have

$H^{\prime}_{con}(1)=-K\sum_{\mathsf{i}}\mathsf{x}_{\mathsf{i}}-\lambda\sum_{\mathsf{p}}\mathcal{Z}_{\mathsf{p}}-\lambda\sum_{\mathsf{p}}\prod_{\mathsf{l}\in\mathsf{p}}Z_{\mathsf{l}}.$
(16)

This is nothing but the usual $\mathbb{Z}_{2}$ gauge theory in its deconfined
phase, stacked with a trivial paramagnet.

Note that the gauge field variables in the condensed phase are precisely the
dipole creation operators. This can be argued intuitively: in the condensed
phase, single fractons can only move with the help of string operators built
from dipoles, and in the gauge theory description these strings become the
Wilson lines which render isolated fractons gauge invariant.

###  III.2 Condensing a subset of fracton dipoles

Identifying the nature of the phases obtained when condensing only a subset of
the fracton dipoles can be done similarly. We give the details here in the
case $N=2$.

####  III.2.1 Condensing a single orientation of dipoles

Suppose first that we condense dipoles with only a single orientation of
dipole moment, which given our conventions is most conveniently chosen to be
${\mathbf{\hat{z}}}$. This can be done simply by using the Hamiltonian (9),
but only including $X_{\mathsf{l}}$, $Z_{\mathsf{l}}$ dipole fields, and the
corresponding hopping terms, for links parallel to ${\mathbf{\hat{z}}}$. We
thus consider the Hamiltonian

$\displaystyle H_{con;z}$
$\displaystyle=-K\sum_{\mathsf{i}}\mathsf{x}_{\mathsf{i}}-g\sum_{\mathsf{c}}(A^{x}_{\mathsf{c}}+A^{y}_{\mathsf{c}})$
(17)
$\displaystyle-\lambda\sum_{\mathsf{p}\perp{\mathbf{\hat{z}}}}Z_{\mathsf{l}_{1}}\mathcal{Z}_{\mathsf{p}}Z_{\mathsf{l}_{2}}\text{,}$

where we have left out $A^{z}_{\mathsf{c}}$ in the second term on the grounds
that it can be written as a product of dipole kinetic energy terms, and where
the last term involves a sum over plaquettes with normal vectors perpendicular
to ${\mathbf{\hat{z}}}$. Note that the pairs $\mathsf{l}_{1},\mathsf{l}_{2}$
are parallel to ${\mathbf{\hat{z}}}$ for both orientations of plaquettes with
$\mathsf{p}\perp{\mathbf{\hat{z}}}$; see figure 7. The Gauss’ law constraint
is an appropriate anisotropic version of (8), viz.

$\mathsf{x}_{\mathsf{i}}\prod_{\mathsf{l}{\,\parallel\,}{\mathbf{\hat{z}}}\ni\mathsf{i}}X_{\mathsf{l}}=\prod_{\mathsf{p}\ni\mathsf{i}}\mathcal{X}_{\mathsf{p}},$
(18)

where the product over $\mathsf{l}$ on the left-hand side is over the two
$z$-axis links $\mathsf{l}$ that touch the vertex $\mathsf{i}$.

We now perform a unitary transformation with the operator
$\mathcal{U}_{z}=\prod_{\mathsf{p}\perp{\mathbf{\hat{z}}}}U_{\mathsf{p}}$,
which is the same as the operator $\mathcal{U}$ defined above, except for the
omission of plaquettes with unit normal along ${\mathbf{\hat{z}}}$. After this
conjugation, we have

$\displaystyle H^{\prime}_{con;z}$
$\displaystyle=\mathcal{U}_{z}H_{con;z}\mathcal{U}_{z}$ (19)
$\displaystyle=-K\sum_{\mathsf{i}}\mathsf{x}_{\mathsf{i}}-g\sum_{\mathsf{c}}(A^{x}_{\mathsf{c}}+A^{y}_{\mathsf{c}})\prod_{\mathsf{l}{\,\parallel\,}{\mathbf{\hat{z}}}\in\mathsf{c}}Z_{\mathsf{l}}-\lambda\sum_{\mathsf{p}\perp{\mathbf{\hat{z}}}}\mathcal{Z}_{\mathsf{p}}\text{.}$

Importantly, the second term is modified by the transformation, which adds a
product of $Z_{\mathsf{l}}$ over the four $z$-axis links $\mathsf{l}$
contained in the cube $\mathsf{c}$. This can clearly be deformed to

$\displaystyle H^{\prime\prime}_{con;z}$
$\displaystyle=-K\sum_{\mathsf{i}}\mathsf{x}_{\mathsf{i}}-g\sum_{\mathsf{c}}[\prod_{\mathsf{p}{\,\parallel\,}{\mathbf{\hat{z}}}\in\mathsf{c}}\mathcal{Z}_{\mathsf{p}}][\prod_{\mathsf{l}{\,\parallel\,}{\mathbf{\hat{z}}}\in\mathsf{c}}Z_{\mathsf{l}}]-\lambda\sum_{\mathsf{p}\perp{\mathbf{\hat{z}}}}\mathcal{Z}_{\mathsf{p}}\text{,}$
(20)

where the first product in the second term is over the two plaquettes on the
boundary of the cube $c$ with normal vectors along ${\mathbf{\hat{z}}}$. Under
conjugation by $\mathcal{U}_{z}$, the Gauss’ law transforms to

$\mathsf{x}_{\mathsf{i}}\prod_{\mathsf{l}{\,\parallel\,}{\mathbf{\hat{z}}}\ni\mathsf{i}}X_{\mathsf{l}}=\prod_{\mathsf{p}{\,\parallel\,}{\mathbf{\hat{z}}}\ni\mathsf{i}}\mathcal{X}_{\mathsf{p}}\text{.}$
(21)

Here the product on the right-hand side is over the four plaquettes with
normals along ${\mathbf{\hat{z}}}$ that meet $\mathsf{i}$ at a corner.

The Hamiltonian (20) together with the above Gauss constraint is nothing but
the “anisotropic model with lineons and planons” discussed in Sec. 5.7 of Ref.
Shirley et al., 2019 (mapped to a generalized gauge theory in the usual way),
with the two types of lineon excitations identified in Ref. Shirley et al.,
2019 corresponding here to single fractons and $\mathfrak{e}^{z}$ lineons.

####  III.2.2 Condensing two orientations of dipoles

Consider now condensing dipoles with dipole moments oriented along
${\mathbf{\hat{z}}}$ and ${\mathbf{\hat{x}}}$. By exchanging dipoles with the
condensate, the fractons in the condensed phase will be able to move along
both the ${\mathbf{\hat{z}}}$ and ${\mathbf{\hat{x}}}$ directions, and hence
the natural expectation is to identify the condensed phase with a decoupled
stack of $\mathbb{Z}_{2}$ gauge theories.

This is indeed what happens. We include $X_{\mathsf{l}}$, $Z_{\mathsf{l}}$ for
$\mathsf{l}$ parallel to ${\mathbf{\hat{z}}}$ and ${\mathbf{\hat{x}}}$, but
not ${\mathbf{\hat{y}}}$. The Gauss’ law is thus

$\mathsf{x}_{\mathsf{i}}\prod_{\mathsf{l}{\,\parallel\,}{\mathbf{\hat{z}}},{\mathbf{\hat{x}}}\ni\mathsf{i}}X_{\mathsf{l}}=\prod_{\mathsf{p}\ni\mathsf{i}}\mathcal{X}_{\mathsf{p}},$
(22)

where the product on the left-hand side is over the $x$ and $z$ links touching
$\mathsf{i}$. For the Hamiltonian we take (9) and include only kinetic terms
for dipoles oriented along the ${\mathbf{\hat{z}}},{\mathbf{\hat{x}}}$
directions. We have, referring to figure 7,

$\displaystyle H_{con;zx}$
$\displaystyle=-K\sum_{\mathsf{i}}\mathsf{x}_{\mathsf{i}}-\lambda\sum_{\mathsf{p}{\,\parallel\,}{\mathbf{\hat{x}}},{\mathbf{\hat{z}}}}Z_{\mathsf{l}_{1}}\mathcal{Z}_{\mathsf{p}}Z_{\mathsf{l}_{2}}$
(23)
$\displaystyle\qquad-\lambda\sum_{\mathsf{p}{\,\parallel\,}{\mathbf{\hat{y}}}}\mathcal{Z}_{\mathsf{p}}(Z_{\mathsf{l}_{1}}Z_{\mathsf{l}_{2}}+Z_{\mathsf{l}_{3}}Z_{\mathsf{l}_{4}}).$

We have again set $h=g=0$, noting that the $A^{a}_{\mathsf{c}}$ terms can be
obtained as products of the dipole kinetic energy terms that are present.
Conjugating with $\mathcal{U}$ as in the analysis of the case where all
dipoles condense, we obtain

$\mathcal{U}H_{con;zx}\mathcal{U}=-K\sum_{\mathsf{i}}\mathsf{x}_{\mathsf{i}}-\lambda\sum_{\mathsf{p}}\mathcal{Z}_{\mathsf{p}}-\lambda\sum_{\mathsf{p}{\,\parallel\,}{\mathbf{\hat{y}}}}\mathcal{Z}_{\mathsf{p}}\prod_{\mathsf{l}\in\mathsf{p}}Z_{\mathsf{l}},$
(24)

while the Gauss law becomes

$\mathsf{x}_{\mathsf{i}}=\prod_{\mathsf{l}{\,\parallel\,}{\mathbf{\hat{z}}},{\mathbf{\hat{x}}}\ni\mathsf{i}}X_{\mathsf{l}}\text{.}$
(25)

After eliminating the now-trivial $\mathcal{Z}_{\mathsf{p}}$ degrees of
freedom, this indeed describes a set of decoupled deconfined $\mathbb{Z}_{2}$
gauge theories, stacked along the ${\mathbf{\hat{y}}}$ direction.

###  III.3 Phase transitions

We now turn to studying the critical points that occur between the
$\mathbb{Z}_{N}$ X-cube phase and the condensed phases described above. We
will first focus on what happens when we only condense fracton dipoles with a
single orientation of dipole moment, which we take to be along
${\mathbf{\hat{z}}}$.

To analyze the critical points we can ignore the fracton matter spins
$\mathsf{x},\mathsf{z}$, as single fractons remain gapped across the
transition. In addition, excitations of the generalized gauge field are gapped
across the transition, so we may also ignore fluctuations of the generalized
gauge field ${\mathcal{Z}}_{\mathsf{p}}$. Indeed, we may set
${\mathcal{Z}}_{\mathsf{p}}=1$ and focus on an effective Hamiltonian for the
dipole matter alone, provided that we remember to only consider correlation
functions of gauge-invariant operators. Since for now we are only interested
in condensing $\mathfrak{d}^{z}_{\mathfrak{m}}$ dipoles, we only need to
retain the corresponding kinetic terms. Thus we can consider the Hamiltonian

$H_{stack}=-h\sum_{\mathsf{l}{\,\parallel\,}{\mathbf{\hat{z}}}}X_{\mathsf{l}}-\lambda\sum_{\mathsf{p}\perp{\mathbf{\hat{z}}}}Z_{\mathsf{l}_{1}}^{\dagger}Z_{\mathsf{l}_{2}}+h.c,$
(26)

where the $\mathsf{l}_{1},\mathsf{l}_{2}$ links are parallel to
${\mathbf{\hat{z}}}$ (see figure 7), and where we are again working with
$\mathbb{Z}_{N}$ spins. Since only the spins on the links parallel to
${\mathbf{\hat{z}}}$ enter into the above Hamiltonian, we see that $H_{stack}$
is simply a decoupled stack of 2d $\mathbb{Z}_{N}$ clock models, with one
clock model at each $z$-coordinate of the lattice.

In this effective matter theory, the $\mathbb{Z}_{N}$ planar conservation laws
governing mobility of lineons (see Sec. II) lead to a subsystem symmetry. In
particular, the number of dipoles in each $xy$ plane is separately conserved
modulo $N$.666Strictly speaking, this conservation only holds for local terms;
non-local processes that add a dipole to every plane are allowed, because each
dipole carries two opposite-sign $\mathbb{Z}_{N}$ charges in neighboring
$xy$-planes. A more general statement is the the difference in number of
dipoles between any two planes is conserved. From a formal perspective, this
symmetry arises from the fact that any local term in the matter theory can be
consistently coupled to the generalized gauge field if and only if it is
symmetry-invariant.

The critical point of the Hamiltonian (26) is then the same as a decoupled
stack of critical 2d $\mathbb{Z}_{N}$ clock model layers, with the added
restriction – coming from gauge invariance – of the effective subsystem
symmetry. The allowed operators are precisely those that are invariant under
independent $\mathbb{Z}_{N}$ transformations on each layer.

####  III.3.1 $N=2$

What happens at the critical point depends on the value of $N$. First,
consider the case of $N=2$. Here the decoupled critical point is built from
Ising∗ CFTs on each layer, where as usual the $*$ denotes the restriction to
the $\mathbb{Z}_{2}$-neutral part of the spectrum. The most relevant way to
couple different layers together is through their energy operators, via a
perturbation to the fixed-point action of the form

$\delta S=\sum_{\alpha\neq\beta}g_{\alpha,\beta}\int
d^{2}x\,d\tau\,\varepsilon_{\alpha}\varepsilon_{\beta},$ (27)

where $\varepsilon_{\alpha}$ is the energy operator on layer $\alpha$. Since
the scaling dimension of the energy operator in the 2d Ising model is
$\Delta_{\varepsilon}\approx 1.41<3/2$,Poland et al. (2019) these couplings
are slightly relevant at the decoupled fixed point. In appendix B we argue
that there is in fact no stable fixed point to this RG flow, with the result
being an instability towards a first-order transition.

####  III.3.2 $N=3$

When $N=3$, we have a 3-state Potts model on each layer. The transition in the
3-state Potts model in three space-time dimensions is well-known to be first
order, and hence we do not obtain any stable fixed points in this case.

####  III.3.3 $N=4$

The case of $N=4$ can be understood by noting that the $\mathbb{Z}_{4}$ clock
model can be mapped to a pair of two decoupled Ising models, written in terms
of two $\mathbb{Z}_{2}$ spins $Z_{1},Z_{2}$ (see e.g. Ortiz et al., 2012).777
To see this, we write the $\mathbb{Z}_{4}$ clock matrix $Z$ in terms of two
$\mathbb{Z}_{2}$ Pauli matrices
$Z_{1}\equiv\sigma^{z}\otimes\mathbf{1},Z_{2}\equiv\mathbf{1}\otimes\sigma^{z}$
as  $Z=\frac{1}{\sqrt{2}}(e^{i\pi/4}Z_{1}+e^{-i\pi/4}Z_{2}).$ (28) The real
matrix $X+X^{\dagger}$ also has a simple representation in terms of
$X_{1}\equiv\sigma^{x}\otimes\mathbf{1}$ and
$X_{2}=\mathbf{1}\otimes\sigma^{x}$, with  $X+X^{\dagger}=X_{1}+X_{2}.$ (29)
It is then easy to check that in terms of the $\mathbb{Z}_{2}$ variables, the
$\mathbb{Z}_{4}$ clock model Hamiltonian splits as $H_{1}+H_{2}$, where
$H_{1}$ ($H_{2}$) is an Ising chain Hamiltonian for the $Z_{1}$ ($Z_{2}$)
variables.  The (gauged) $\mathbb{Z}_{4}$ symmetry acts in this representation
as

$\mathbb{Z}_{4}\,:\,Z_{1}\rightarrow Z_{2},\qquad Z_{2}\rightarrow
Z_{1}X_{1},$ (30)

and as such at the critical point it exchanges the two energy operators
$\varepsilon_{1}$ and $\varepsilon_{2}$. The operator
$\varepsilon_{1}+\varepsilon_{2}$ is therefore gauge-invariant, and since the
Ising models are decoupled $\varepsilon_{1}+\varepsilon_{2}$ has the same
dimension as the energy operator in a single Ising model. Therefore according
to the discussion of the $N=2$ case above, the decoupled fixed point is again
unstable.

####  III.3.4 $N>4$

When $N>4$, we describe the critical point on each layer with the action

$\displaystyle S$ $\displaystyle=\int
d\tau\,d^{2}x\,\Big{(}|\partial\psi|^{2}+t|\psi|^{2}+\frac{u}{4}|\psi|^{4}$
(31)
$\displaystyle\qquad+\frac{g\Lambda^{4-N}}{N!}(\psi^{N}+(\psi^{*})^{N})+\cdots\Big{)},$

with $\psi$ a complex scalar field. It is known from Monte Carlo simulations
that the $\mathbb{Z}_{N}$ anisotropy term is irrelevant as long as
$N>4$,Scholten and Irakliotis (1993); Hove and Sudbo (2003) so that the phase
transitions on each layer are in the universality class of the 2+1D XY∗ model.
The most relevant gauge-invariant coupling between the layers is the energy-
energy coupling. The energy operator in this case is known to have scaling
dimensionPoland et al. (2019)

$\Delta_{\varepsilon}\approx 1.51>3/2,$ (32)

and so the decoupled fixed point consisting of a decoupled stack of 2+1D XY∗
models is (barely) stable.

####  III.3.5 No particle-hole symmetry

For $N>2$, the discussion above strictly speaking pertains only to models
which possess a $\mathbb{Z}_{2}$ particle-hole symmetry $C$, which acts on the
spins of the $\mathbb{Z}_{N}$ clock models (as well as on the original spins
of the XC model) by conjugation with the $N\times N$ matrix $C$, where

$C=\begin{pmatrix}1&&&\\\ &&&1\\\ &&\scalebox{-1.0}[1.0]{$\ddots$}&\\\
&1&&\end{pmatrix},\qquad C\mathsf{z}C=\mathsf{z}^{\dagger},\quad
C\mathsf{x}C=\mathsf{x}^{\dagger}.$ (33)

In terms of the coarse-grained field $\psi$ employed above in (31), this sends
$C:\psi\rightarrow\psi^{*}$.

This symmetry can be broken for $N>2$ by letting the parameters $h,\lambda$ in
(26) be complex. We are then prompted to consider an action on each layer of
the form

$\displaystyle S$ $\displaystyle=\int
d\tau\,d^{2}x\,\Big{(}\psi^{*}\partial_{\tau}\psi+\frac{1}{2m}|\nabla\psi|^{2}+t|\psi|^{2}+\frac{u}{4}|\psi|^{4}$
(34)
$\displaystyle\qquad+\frac{g\Lambda^{4-N}}{N!}(\psi^{N}+(\psi^{*})^{N})+\cdots\Big{)},$

where the $C$-breaking first term, which leads to non-relativistic $z=2$
scaling, is included instead of $|\partial_{\tau}\psi|^{2}$ on the grounds of
it being more relevant. Since 2+1D is the upper critical dimension for the
above model, the RG flow can be accurately computed through a perturbative
analysis of the above action (34).

The non-relativistic nature of the theory means that the 1-loop computation of
the $\beta$ functions is exact,888For example, in momentum-shell
regularization with a cutoff $\Lambda$, the only diagrams depending on
$\ln(\Lambda)$ are those built from concatenations of one-loop bubbles,
meaning that the full beta functions are determined by the one-loop terms.
giving

$\beta_{u}=-u^{2}/2,\qquad\beta_{g}=(4-N)g-ug/2.$ (35)

The anisotropy is therefore relevant at $N=3$, marginally irrelevant at $N=4$,
and irrelevant for $N>4$. When $N=3$ Ginzburg-Landau theory predicts a first-
order transition due to the shape of the potential for $\psi$, and so this
case can be ignored. For $N\geq 4$ however, we obtain a nontrivial critical
theory describable by the non-relativistic XY∗ model.

We therefore need to understand the stability of a decoupled stack of non-
relativistic XY∗ models.999This story is very similar to the analysis of
stacks of Fermi liquid $\leftrightarrow$ orthogonal metal phase transitions
studied in Ref. Zou and Senthil (2016) The most relevant gauge-invariant
couplings are again those which couple the energy operators on different
layers; in the present notation they read

$\delta S=\sum_{\alpha,\beta}g_{\alpha,\beta}\int
d\tau\,d^{2}x\,|\psi_{\alpha}|^{2}|\psi_{\beta}|^{2}.$ (36)

The beta function for $g_{\alpha,\beta}$ is similarly exactly computable from
the 1-loop term, which gives

$\beta_{g_{\alpha,\beta}}=-\frac{1}{2}(g_{\alpha,\beta})^{2}.$ (37)

Due to the non-relativistic nature of the theory, none of the
$g_{\alpha,\beta}$s mix with each other. As a result, the decoupled fixed
point is stable, provided that all of the $g_{\alpha,\beta}$s are
positive.101010Note that getting a stable fixed point does not require fine-
tuning: we just need to that $g_{\alpha,\beta}>0$ for all $\alpha,\beta$; we
do not require that each of the $g_{\alpha,\beta}$ be tuned to any one
particular value.

Before moving on, we add a brief aside about a subtlety in the above
discussion of stability. In the argument above, we have restricted our
attention only to the $g_{\alpha,\beta}$ couplings, which are marginal. In
principle however, the irrelevant operators we have neglected may drive some
of the $g_{\alpha,\beta}$ negative during the initial stages of the RG flow.
If this occurs generically, the stability of the proposed fixed point would be
called into question.

This worry is not a serious issue in the present context, however. First, note
that due to the non-relativistic scaling, the only irrelevant operators which
can renormalize $g_{\alpha,\beta}$ are those containing
$\psi^{*}_{\alpha}\psi_{\alpha}\psi^{*}_{\beta}\psi_{\beta}$ and no other
field operators, the least irrelevant of which is
$(\psi_{\alpha}^{*}\nabla_{i}\psi_{\alpha})(\psi_{\beta}^{*}\nabla^{i}\psi_{\beta})$,
which has an RG eigenvalue of $-2$. Unless the bare value of this term is
significantly larger than that of $g_{\alpha,\beta}$, $g_{\alpha,\beta}$ will
never be driven negative during the flow. Therefore if we make the
(reasonable) assumption that the bare values of these irrelevant operators are
not significantly larger than the bare values of the appropriate
$g_{\alpha,\beta}$, the above fixed point is indeed stable.111111It is
instructive to compare the story here with the Kohn-Luttinger instability that
occurs in Fermi liquids.Shankar (1994) In that scenario, irrelevant operators
lead to instabilities by renormalizing certain marginal couplings. This occurs
however only due to the fact that the bare values of certain marginal terms
are naturally exponentially smaller than the bare values of the irrelevant
operators which renormalize them. In the present setting, by contrast, the
bare values of the irrelevant terms and of the marginal terms they renormalize
both scale similarly with the distance between the $\alpha$ and $\beta$
layers.

####  III.3.6 Condensing multiple species of dipoles

So far we have only been concerned with condensing a single species of dipole,
but the same analysis can be applied when simultaneously condensing multiple
species. The starting point for the critical theory is simply multiple
decoupled stacks of critical $\mathbb{Z}_{N}$ clock models, with one stack
oriented along the dipole moment vector of each of the condensing dipoles. The
effective matter theory again has a subsystem symmetry, where (for local
terms) the number of $z$-dipoles is conserved modulo $N$ in each $xy$-plane,
and similarly for $x$\- and $y$-dipoles. Dipoles with different orientations
of dipole moment are separately conserved, because they carry $\mathbb{Z}_{N}$
planar charges in distinct planes. The discussion of stability of the
decoupled critical point is not modified, because the couplings between
perpendicular planes are always less relevant than the couplings between
parallel planes.

##  IV Lineon dipole condensation

In this section we discuss what happens when lineon dipoles are condensed. The
analysis is similar in many respects to the case of fracton dipoles treated in
the previous section, although the phenomenology is slightly richer due to
fact that there are three different types of lineons, as opposed to only a
single type of fracton. As in the previous section, we will proceed by first
using a generalized gauge theory to identify the nature of the condensed
phases, and then later turn to discussing the nature of the phase transitions.

###  IV.1 Condensed phases and dual gauge theories

Figure 8:  An illustration of a string operator moving an $\mathfrak{e}^{x}$
excitation in the ${\mathbf{\hat{z}}}$ direction in the presence of a
$\mathfrak{d}^{z}_{\mathfrak{e}}$ condensate. The extra excitations created by
the middle part of the string may all be absorbed into the condensate.

Before discussing concrete Hamiltonians, let us assess our expectations for
the condensed phases. First consider condensing a single species of dipole,
say $\mathfrak{d}^{z}_{\mathfrak{e}}$. In the electric sector,
$\mathfrak{e}^{x}$ and $\mathfrak{e}^{y}$ lineons now become capable of moving
along the ${\mathbf{\hat{z}}}$ direction, as they can do so by absorbing
dipoles from the condensate. An example of a string operator implementing this
type of process is shown in figure 8, where an $\mathfrak{e}^{x}$ lineon moves
in the ${\mathbf{\hat{z}}}$ direction by absorbing
$\mathfrak{d}^{z}_{\mathfrak{e}}$ dipoles from the condensate. The
$\mathfrak{e}^{x}$ and $\mathfrak{e}^{y}$ lineons still remain distinct
excitations however, as there continues to be no local process which converts
an isolated $\mathfrak{e}^{x}$ into an isolated $\mathfrak{e}^{y}$. The
electric sector thus consists of two types of $\mathbb{Z}_{N}$ charges, each
of which are free to move in planes normal to the ${\mathbf{\hat{x}}}$ and
${\mathbf{\hat{y}}}$ directions.

Figure 9: Processes which give tension to fracton membranes in the presence of
a $\mathfrak{d}^{z}_{\mathfrak{e}}$ condensate. In both panels, purple links
denote the support of fracton membrane operators, while orange links denote
the support of condensed lineon dipole strings. a) a process giving line
tension to a membrane in an $xz$ plane. b) a corner-turning process giving
surface tension to a membrane in an $xy$ plane.

In the magnetic sector, the condensate gives a tension to the membrane
operators creating fractons at their corners in a rather subtle way. For such
membrane operators lying in a $xz$-plane, the string operators of the
condensed $\mathfrak{d}^{z}_{\mathfrak{e}}$ dipoles only fail to commute with
the membrane near its edges running along the $x$-direction; an example of
such a processes is shown in the left panel of Fig. 9. This means that the
$x$-edges of the membrane acquire a line tension. Similarly, the $y$-edges of
membranes lying in a $yz$-plane also get a line tension from the condensate.
On the other hand, membranes lying in $xy$-planes get a surface tension,
because the string operators of condensate dipoles lying in the same
$xy$-plane anticommute with the membrane operator, if the dipole turns a
corner anywhere within the area of the membrane; an example is shown in the
right panel of Fig. 9. These results imply that isolated fractons, as well as
$\mathfrak{d}^{z}_{\mathfrak{m}}$ dipoles are confined. However the
$\mathfrak{d}^{x,y}_{\mathfrak{m}}$ dipoles remain deconfined.

The deconfined excitations of the theory are thus the same as that of two
stacks of deconfined 2d $\mathbb{Z}_{N}$ gauge theories, with one stack
oriented along ${\mathbf{\hat{y}}}$ and another along ${\mathbf{\hat{x}}}$. We
will see later that the condensed theory is indeed given by two
interpenetrating stacks of 2d $\mathbb{Z}_{N}$ gauge theories.

Now consider condensing a second species of dipole, say
$\mathfrak{d}^{y}_{\mathfrak{e}}$. Since $\mathfrak{d}^{y}_{\mathfrak{e}}$ is
a bound state of the gauge charges in two adjacent layers in the gauge theory
stack oriented along ${\mathbf{\hat{y}}}$, condensing
$\mathfrak{d}^{y}_{\mathfrak{e}}$ will turn this stack into a single 3d gauge
theory, while leaving the stack oriented along ${\mathbf{\hat{x}}}$ unchanged.
We therefore obtain a 3d $\mathbb{Z}_{N}$ gauge theory (with gauge charge
$\mathfrak{e}^{x}$) and a stack of 2d $\mathbb{Z}_{N}$ gauge theories oriented
along the ${\mathbf{\hat{x}}}$ direction (with gauge charges given by
$\mathfrak{e}^{y}$ on each layer).

Finally, consider the phase obtained when all the
$\mathfrak{d}^{a}_{\mathfrak{e}}$s are condensed. Each $\mathfrak{e}^{a}$ is
now free to move in all three directions, and since there is still no local
process turning an $\mathfrak{e}^{x}$ into an $\mathfrak{e}^{y}$, a natural
assumption is that the condensed phase is a deconfined 3d $\mathbb{Z}_{N}^{2}$
gauge theory.

The above logic gives plausible identifications for the condensed phases in
each case, but is not completely rigorous. For example, while the electric
sector in the phase where all $\mathfrak{d}^{a}_{\mathfrak{e}}$ condense is
the same as that of $\mathbb{Z}_{N}^{2}$ gauge theory, it is not obvious how
the magnetic degrees of freedom in the XC model end up organizing themselves
into a $\mathbb{Z}_{N}^{2}$’s worth of loop excitations.

To conclusively demonstrate that these expectations are borne out, we will
employ a generalized gauge theory construction similar to that used when
discussing the condensation of fracton dipoles.

We will first focus on the case where all orientations of dipoles condense. We
will furthermore specialize to the $\mathbb{Z}_{2}$ case for simplicity of
notation; the generalization to $\mathbb{Z}_{N}$ is done using the same
methods as employed in Appendix A. According to the discussion above, we
expect that the resulting phase will be a deconfined $\mathbb{Z}_{2}^{2}$
gauge theory. To show this, it is helpful to realize the XC model as a
particular limit of a Hamiltonian defined in a larger Hilbert space. Since we
are aiming for a $\mathbb{Z}_{2}^{2}$ gauge theory, we will want to have two
degrees of freedom on each link. In fact we will find it more convenient to
work with a model containing three qubits on each link, together with a single
constraint reducing the total number of degrees of freedom on each link to
two.

Figure 10:  Representations of the gauge field terms appearing in the
generalized Gauss laws. The operators $A^{ab}$ are defined as products of
pairs of $X^{c}X^{d}$ operators over each of the black links, with $c,d$
indicated by the pairs of letters.

We thus start with a Hamiltonian defined on a Hilbert space with three qubits
on each link of the direct lattice, with Pauli operators $X_{l}^{a},Z_{l}^{a}$
($a=x,y,z$), which are subject to the constraint

$Z^{x}_{l}Z^{y}_{l}Z^{z}_{l}=1$ (38)

for each link $l$. Within this Hilbert space, we can write a Hamiltonian that
reduces to $H_{XC}$ as

$H_{XC}^{\prime}=-g\sum_{i}\sum_{a\neq b}A_{i}^{ab}-K\sum_{c}\prod_{l\in
c}Z^{l}_{l}-h\sum_{l}X^{a}_{l}X^{b}_{l},$ (39)

where $Z^{l}_{l}$ means $Z^{a}_{l}$ with ${\mathbf{\hat{a}}}{\,\parallel\,}l$,
and where in the last term $a\neq b$ with
${\mathbf{\hat{a}}},{\mathbf{\hat{b}}}\perp l$. The operators $A^{ab}_{i}$,
$a\neq b$ are defined as

$\displaystyle A^{ab}_{i}=\prod_{l\perp{\mathbf{\hat{a}}}\ni
i}X^{b}_{l}\prod_{l^{\prime}\perp{\mathbf{\hat{b}}}\ni
i}X^{a}_{l^{\prime}}\prod_{l^{\prime\prime}\perp{\mathbf{\hat{c}}}\ni
i}X^{c}_{l^{\prime\prime}},$ (40)

where $c$ is neither $a$ nor $b$, and where for instance the first product is
over the four links touching $i$ that are perpendicular to
${\mathbf{\hat{a}}}$. Note that $A^{ab}_{i}=A^{ba}_{i}$. An illustration of
the three different types of $A^{ab}$ operators is given in figure 10.
$H_{XC}^{\prime}$ preserves the constraint (38), since each term contains an
even number of $X^{a}_{l}$ operators on each link. Moreover, all the terms in
$H_{XC}^{\prime}$ commute with one another.

To see why $H_{XC}^{\prime}$ is equivalent to $H_{XC}$, we work in the ground
state subspace of the $h$-term, where we have the additional local constraint
$X^{a}_{l}X^{b}_{l}=1$ on each link. (Again, $a\neq b$ and
${\mathbf{\hat{a}}},{\mathbf{\hat{b}}}\perp l$.) We thus have a single
effective qubit on each link, with Pauli operators

$X^{{\rm eff}}_{l}=X^{c}_{l}X^{a}_{l}=X^{c}_{l}X^{b}_{l}$ (41)

and

$Z^{{\rm eff}}_{l}=Z^{c}_{l}=Z^{a}_{l}Z^{b}_{l}\text{,}$ (42)

where ${\mathbf{\hat{c}}}{\,\parallel\,}l$. The remaining $g$\- and $K$-terms
of $H_{XC}^{\prime}$ are easily expressed in terms of $X^{{\rm eff}}_{l}$ and
$Z^{{\rm eff}}_{l}$, and reduce to $H_{XC}$.

As in the case of fracton dipole condensation, in order to condense lineon
dipoles while keeping single lineons gapped, it is convenient first to map
$H^{\prime}_{XC}$ to a generalized gauge theory. The gauge theory is
constructed by placing three qubits on each vertex and on each link of the
lattice. We first discuss the site variables, whose Pauli operators are
denoted $\mathsf{x}^{a}_{i},\mathsf{z}^{a}_{i}$ on each vertex $i$, and
$a=x,y,z$. We will see that $\mathsf{z}^{a}_{i}$ creates an $\mathfrak{e}^{a}$
lineon at the vertex $i$. In accordance with the lineon fusion rules, we
impose the constraint

$\mathsf{z}^{x}_{i}\mathsf{z}^{y}_{i}\mathsf{z}^{z}_{i}=1$ (43)

at each site.

On each link, the Hilbert space is the same as in the construction of the
ungauged model $H^{\prime}_{XC}$, except that we denote the Pauli operators
for the three qubits by $\widetilde{X}_{l}^{a},\widetilde{Z}^{a}_{l}$. We
impose the same constraint on each link, namely
$\widetilde{Z}^{x}_{l}\widetilde{Z}^{y}_{l}\widetilde{Z}^{z}_{l}=1$. We define
operators $\widetilde{A}^{ab}_{i}$ in terms of the $\widetilde{X}_{l}$ by the
same formula as (40).

To reduce down to the original number of degrees of freedom, we then impose,
for each $a\neq b$, the Gauss law constraints

$\displaystyle\mathsf{x}_{i}^{a}\mathsf{x}_{i}^{b}=\widetilde{A}_{i}^{ab}.$
(44)

Note that this is compatible with the constraint (43). Moreover, only two of
the three equations in (44) are independent, because taking the product of all
three equations gives the triviality $1=1$.

The mapping is established by the following dictionary between the operators
in $H^{\prime}_{XC}$ and those in the generalized gauge theory:

$\displaystyle X^{a}_{l}X^{b}_{l}$ $\displaystyle\mapsto$
$\displaystyle\widetilde{X}^{a}_{l}\widetilde{X}^{b}_{l}\qquad(\text{any pair
}a\neq b)$ (45) $\displaystyle Z^{a}_{{\langle ij\rangle}}$
$\displaystyle\mapsto$
$\displaystyle\mathsf{z}^{a}_{i}\widetilde{Z}^{a}_{{\langle
ij\rangle}}\mathsf{z}^{a}_{j}\qquad({\mathbf{\hat{a}}}{\,\parallel\,}{\langle
ij\rangle})$ (46) $\displaystyle Z^{b}_{{\langle ij\rangle}}$
$\displaystyle\mapsto$
$\displaystyle\mathsf{z}^{c}_{i}\widetilde{Z}^{b}_{{\langle
ij\rangle}}\mathsf{z}^{c}_{j}\qquad({\mathbf{\hat{b}}},{\mathbf{\hat{c}}}\perp{\langle
ij\rangle},b\neq c)\text{.}$ (47)

Under this mapping, the modified X-cube Hamiltonian $H^{\prime}_{XC}$ becomes

$\widetilde{H}^{\prime}_{XC}=-g\sum_{i}\sum_{a\neq
b}\mathsf{x}^{a}_{i}\mathsf{x}^{b}_{i}-K\sum_{c}\prod_{l\in
c}\widetilde{Z}^{l}_{l}-h\sum_{l}\widetilde{X}^{a}_{l}\widetilde{X}^{b}_{l}\text{.}$
(48)

The operator dictionary tells us that, on a link ${\langle
ij\rangle}{\,\parallel\,}{\mathbf{\hat{a}}}$, the non-gauge-invariant operator
$\widetilde{Z}^{b}_{\langle ij\rangle}$ creates an excitation in the same
superselection sector as the $\mathfrak{d}^{z}_{\mathfrak{e}}$ lineon dipole
$\mathsf{z}^{c}_{i}\mathsf{z}^{c}_{j}$, where
${\mathbf{\hat{b}}},{\mathbf{\hat{c}}}\perp{\langle ij\rangle}$ and $b\neq c$.
This is so because the gauge-invariant local operator
$\mathsf{z}^{c}_{i}\widetilde{Z}^{b}_{{\langle ij\rangle}}\mathsf{z}^{c}_{j}$
necessarily creates a trivial, locally createable excitation. We also observe
that acting with $\widetilde{Z}^{b}_{\langle ij\rangle}$ violates the $h$-term
in $\widetilde{H}^{\prime}_{XC}$, and no other terms, so that $h$ can be
understood as setting the energy gap for lineon dipoles. Since $g$ controls
the gap for single lineons, we can control the single-lineon and lineon-dipole
gaps independently in this formulation, as desired. We note that, by the same
reasoning as above, acting with $\widetilde{Z}^{a}_{\langle ij\rangle}$
(${\mathbf{\hat{a}}}{\,\parallel\,}{\langle ij\rangle}$) creates an excitation
in the same sector as a pair of $\mathfrak{e}^{a}$ lineons separated in the
${\mathbf{\hat{a}}}$-direction; however, as discussed in Sec. II, such a pair
is a trivial excitation.

Figure 11: An illustration of the operators $B^{ab}_{p}$. The indices on each
link $l$ denote the corresponding $\widetilde{Z}^{a}_{l}$ operators.

Now we introduce operators that, when added to the Hamiltonian, will play the
role of kinetic energy for the lineon dipoles. Recall that a lineon dipole is
a planon, moving in the plane perpendicular to its dipole moment. On each
plaquette $p$ lying in the $ab$-plane, we define two operators, $B^{ab}_{p}$
and $B^{ba}_{p}$. $B^{ab}_{p}$ hops a $\mathfrak{d}^{a}_{\mathfrak{e}}$ dipole
along ${\mathbf{\hat{b}}}$, where the dipole consists of a pair of
$\mathfrak{e}^{b}$ excitations. Similarly, $B^{ba}_{p}$ hops a
$\mathfrak{d}^{b}_{\mathfrak{e}}$ (consisting of a pair of
$\mathfrak{e}^{a}$’s) along ${\mathbf{\hat{a}}}$. More concretely, we define

$\displaystyle B^{zx}_{p}$ $\displaystyle=$
$\displaystyle\widetilde{Z}^{y}_{l_{1}}\widetilde{Z}^{x}_{l_{3}}\widetilde{Z}^{x}_{l_{4}}\widetilde{Z}^{y}_{l_{2}}$
(49) $\displaystyle B^{xz}_{p}$ $\displaystyle=$
$\displaystyle\widetilde{Z}^{y}_{l_{3}}\widetilde{Z}^{z}_{l_{1}}\widetilde{Z}^{z}_{l_{2}}\widetilde{Z}^{y}_{l_{4}}$
(50) $\displaystyle B^{zy}_{p}$ $\displaystyle=$
$\displaystyle\widetilde{Z}^{x}_{l_{1}}\widetilde{Z}^{y}_{l_{3}}\widetilde{Z}^{y}_{l_{4}}\widetilde{Z}^{x}_{l_{2}}$
(51) $\displaystyle B^{yz}_{p}$ $\displaystyle=$
$\displaystyle\widetilde{Z}^{x}_{l_{3}}\widetilde{Z}^{z}_{l_{1}}\widetilde{Z}^{z}_{l_{2}}\widetilde{Z}^{x}_{l_{4}}$
(52) $\displaystyle B^{xy}_{p}$ $\displaystyle=$
$\displaystyle\widetilde{Z}^{z}_{l_{1}}\widetilde{Z}^{y}_{l_{3}}\widetilde{Z}^{y}_{l_{4}}\widetilde{Z}^{z}_{l_{2}}$
(53) $\displaystyle B^{yx}_{p}$ $\displaystyle=$
$\displaystyle\widetilde{Z}^{z}_{l_{3}}\widetilde{Z}^{x}_{l_{1}}\widetilde{Z}^{x}_{l_{2}}\widetilde{Z}^{z}_{l_{4}}\text{.}$
(54)

These operators are displayed graphically in Fig. 11. The indexing of links in
each plaquette is exactly as illustrated for dual lattice plaquettes in Fig.
7. These operators are easily checked to be gauge-invariant. In each term, the
two factors of $\widetilde{Z}^{a}_{l}$ written on the outside of the product
create/destroy lineon dipoles, while the two factors of
$\widetilde{Z}^{a}_{l}$ written on the inside are needed to make the term
gauge-invariant. We see that for ${\mathbf{\hat{a}}}{\,\parallel\,}l$,
$\widetilde{Z}^{a}_{l}$ plays the role of a gauge field, while for
${\mathbf{\hat{a}}}\perp l$, $\widetilde{Z}^{a}_{l}$ is a matter field
creating a lineon dipole. By taking products for instance of $B^{zx}_{p}$ and
$B^{zy}_{p}$ operators, we can obtain string operators for a
$\mathfrak{d}^{z}_{\mathfrak{e}}$ dipole moving along an arbitrary path in the
$xy$-plane. Whenever the dipole turns a corner, the string operator contains a
factor of $\widetilde{Z}^{z}_{l}$ with $l{\,\parallel\,}{\mathbf{\hat{z}}}$,
corresponding to the two constituent particles converting from
$\mathfrak{e}^{x}$’s to $\mathfrak{e}^{y}$’s (or vice versa) by exchanging an
$\mathfrak{e}^{z}$.

The following solvable Hamiltonian condenses all species of lineon dipoles,
while keeping single lineons gapped:

$H_{con}=-g\sum_{i}\sum_{a\neq
b}\mathsf{x}^{a}_{i}\mathsf{x}^{b}_{i}-\lambda\sum_{a\neq
b}\sum_{p{\,\parallel\,}{\mathbf{\hat{a}}}\times{\mathbf{\hat{b}}}}(B^{ab}_{p}+B^{ba}_{p})\text{.}$
(55)

Here, in the last term, the sumis over plaquettes $p$ lying in the $ab$-plane.
In addition to adding the kinetic energy term to
$\widetilde{H}^{\prime}_{XC}$, we have set $h=0$ to allow the lineon dipoles
to condense and obtain a solvable Hamiltonian. Moreover, we have also set
$K=0$, as the $K$-term can be obtained as a product of the $B^{ab}_{p}$
kinetic energy operators.

It turns out that $H_{con}$ is the $\mathbb{Z}_{2}^{2}$ gauge theory we
expected to find. To see this, we first recall that given the constraint
$\widetilde{Z}^{x}_{l}\widetilde{Z}^{y}_{l}\widetilde{Z}^{z}_{l}=1$, we
effectively have two qubits on each link. We define $Z^{1,2}_{l},X^{1,2}_{l}$
Pauli operators for these two qubits, where the upper index will label the two
$\mathbb{Z}_{2}$ gauge theories. Specifically we choose

$\displaystyle Z^{1}_{l}=\widetilde{Z}^{z}_{l},$ $\displaystyle
Z^{2}_{l}=\widetilde{Z}^{y}_{l}\qquad$
$\displaystyle(l{\,\parallel\,}{\mathbf{\hat{x}}})$ (56) $\displaystyle
Z^{1}_{l}=\widetilde{Z}^{y}_{l},$ $\displaystyle
Z^{2}_{l}=\widetilde{Z}^{x}_{l}\qquad$
$\displaystyle(l{\,\parallel\,}{\mathbf{\hat{y}}})$ (57) $\displaystyle
Z^{1}_{l}=\widetilde{Z}^{x}_{l},$ $\displaystyle
Z^{2}_{l}=\widetilde{Z}^{z}_{l}\qquad$
$\displaystyle(l{\,\parallel\,}{\mathbf{\hat{z}}})\text{.}$ (58)

From this choice it follows that
$X^{1}_{l}=\widetilde{X}^{z}_{l}\widetilde{X}^{x}_{l}$ and
$X^{2}_{l}=\widetilde{X}^{y}_{l}\widetilde{X}^{x}_{l}$ for
$l{\,\parallel\,}{\mathbf{\hat{x}}}$, with similar choices for the two other
orientations of $l$. For the lineon matter fields, we define

$\displaystyle\mathsf{x}^{1}_{i}=\mathsf{x}^{x}_{i}\mathsf{x}^{y}_{i}\qquad$
$\displaystyle\mathsf{x}^{2}_{i}=\mathsf{x}^{x}_{i}\mathsf{x}^{z}_{i}$ (59)
$\displaystyle\mathsf{z}^{1}_{i}=\mathsf{z}^{y}_{i}\qquad$
$\displaystyle\mathsf{z}^{2}_{i}=\mathsf{z}^{z}_{i}\text{.}$ (60)

With these choices the Gauss law equations for $ab=xy$ and $ab=xz$ become

$\displaystyle\widetilde{A}^{xy}_{i}$ $\displaystyle=$ $\displaystyle
A^{1}_{i}=\mathsf{x}^{1}_{i}$ (61) $\displaystyle\widetilde{A}^{xz}_{i}$
$\displaystyle=$ $\displaystyle A^{2}_{i}=\mathsf{x}^{2}_{i}\text{,}$ (62)

where $A^{1,2}_{i}=\prod_{l\ni i}X^{1,2}_{l}$ is a product over the six links
touching $i$. Similarly defining $B^{1,2}_{p}=\prod_{l\in p}Z^{1,2}_{l}$, the
Hamiltonian becomes

$\displaystyle H_{con}$ $\displaystyle=$
$\displaystyle-g\sum_{i}(\mathsf{x}^{1}_{i}+\mathsf{x}^{2}_{i})-\lambda\sum_{p{\,\parallel\,}{\mathbf{\hat{x}}}}(B^{1}_{p}+B^{2}_{p})$
(63) $\displaystyle-$
$\displaystyle\lambda\sum_{p{\,\parallel\,}{\mathbf{\hat{y}}}}(B^{2}_{p}+B^{1}_{p}B^{2}_{p})-\lambda\sum_{p{\,\parallel\,}{\mathbf{\hat{z}}}}(B^{1}_{p}+B^{1}_{p}B^{2}_{p})\text{.}$

This is not quite the standard Hamiltonian for two decoupled $\mathbb{Z}_{2}$
gauge theories, due to the presence of the $B^{1}_{p}B^{2}_{p}$ terms.
However, acting on the ground state we clearly have $B^{1}_{p}=B^{2}_{p}=1$,
so the ground state is the same as that of the standard Hamiltonian. Moreover,
it is trivial to deform this Hamiltonian to the standard one by the obvious
linear interpolation, so $H_{con}$ is indeed in the same phase as the standard
deconfined $\mathbb{Z}_{2}^{2}$ gauge theory.

###  IV.2 Condensing a subset of lineon dipoles

We can similarly consider what happens when only a subset of the lineon
dipoles condense. First consider condensing only a single orientation of
dipole, say with dipole moment along ${\mathbf{\hat{z}}}$. To treat this case,
we modify (55) as follows:

$\displaystyle H_{con;z}$ $\displaystyle=$ $\displaystyle-g\sum_{i}\sum_{a\neq
b}\mathsf{x}^{a}_{i}\mathsf{x}^{b}_{i}-\lambda\sum_{p{\,\parallel\,}{\mathbf{\hat{z}}}\times{\mathbf{\hat{x}}}}B^{zx}_{p}-\lambda\sum_{p{\,\parallel\,}{\mathbf{\hat{z}}}\times{\mathbf{\hat{y}}}}B^{zy}_{p}$
(64) $\displaystyle-$ $\displaystyle
h\sum_{l{\,\parallel\,}{\mathbf{\hat{x}}}}\widetilde{X}^{y}_{l}\widetilde{X}^{z}_{l}-h\sum_{l{\,\parallel\,}{\mathbf{\hat{y}}}}\widetilde{X}^{x}_{l}\widetilde{X}^{z}_{l}\text{.}$

Here we include the kinetic energy terms only for
$\mathfrak{d}^{z}_{\mathfrak{e}}$ dipoles, and retain the $h$-terms on $x$\-
and $y$-links, which keep $\mathfrak{d}^{x,y}_{\mathfrak{e}}$ dipoles gapped.
The $K$-term is not included as it can be obtained as a product of kinetic
energy terms. A similar analysis to that above shows that, in the ground state
subspace of the $h$-term, $H_{con;z}$ is precisely two independent stacks of
2d $\mathbb{Z}_{2}$ gauge theories, stacked along the ${\mathbf{\hat{x}}}$ and
${\mathbf{\hat{y}}}$ directions.

Now consider condensing dipoles with moments along both ${\mathbf{\hat{z}}}$
and ${\mathbf{\hat{x}}}$. In this case we consider

$\displaystyle H_{con;z}$ $\displaystyle=$ $\displaystyle-g\sum_{i}\sum_{a\neq
b}\mathsf{x}^{a}_{i}\mathsf{x}^{b}_{i}-\lambda\sum_{p{\,\parallel\,}{\mathbf{\hat{z}}}\times{\mathbf{\hat{x}}}}(B^{zx}_{p}+B^{xz}_{p})$
$\displaystyle-$
$\displaystyle\lambda\sum_{p{\,\parallel\,}{\mathbf{\hat{z}}}\times{\mathbf{\hat{y}}}}B^{zy}_{p}-\lambda\sum_{p{\,\parallel\,}{\mathbf{\hat{x}}}\times{\mathbf{\hat{y}}}}B^{xy}_{p}-h\sum_{l{\,\parallel\,}{\mathbf{\hat{y}}}}\widetilde{X}^{x}_{l}\widetilde{X}^{z}_{l}\text{.}$

This includes kinetic energy terms for $\mathfrak{d}^{x,z}_{\mathfrak{e}}$
dipoles, while keeping $\mathfrak{d}^{y}_{\mathfrak{e}}$ dipoles gapped.
Again, a similar analysis shows that this model is in the same phase as a 3d
$\mathbb{Z}_{2}$ gauge theory plus a stack along ${\mathbf{\hat{y}}}$ of 2d
$\mathbb{Z}_{2}$ gauge theories, as expected.

The generalization of these constructions to the $\mathbb{Z}_{N}$ case is
straightforward, and can be done in much the same manner as described in
appendix A for the case of fracton dipoles. All the prescriptions for adding
Hermitian conjugates to operators are the same as in that case, and as such we
will omit the details.

###  IV.3 Phase transitions

Because we are again condensing excitations that can move in 2d planes, with
the number of excitations in each plane separately conserved modulo $N$ (for
local processes), the critical points to analyze are the same as those that
appeared in section III.3. The stability analysis of the critical points is
identical to the previous fracton dipole case, and so we will not repeat
ourselves.

##  V Lineon condensation

In this section we will discuss what happens when we condense one or more
species of the lineons $\mathfrak{e}^{x,y,z}$. We will see that when we
condense a single species of lineon the resulting phase is a stack of
deconfined 2d $\mathbb{Z}_{N}$ gauge theories, while if more species are
condensed the resulting phase is trivial. If $N>4$, we find that the system
enters an intermediate gapless phase upon condensation of lineons, described
as an array of strongly-coupled Luttinger liquids. The transition from the XC
phase into the gapless phase is expected to be of Berezinskii-Kosterlitz-
Thouless (BKT) type.

###  V.1 Condensed phases and dual Hamiltonians

Before getting into details let us discuss our expectations for the condensed
phases. Consider first condensing a single type of lineon, for example
$\mathfrak{e}^{z}$. In the electric sector, the fusion rule
$\mathfrak{e}^{x}\mathfrak{e}^{y}\mathfrak{e}^{z}=1$ means that the condensed
phase contains a single type of electric excitation
$\mathfrak{e}^{x}\sim(\mathfrak{e}^{y})^{\dagger}$, which is able to move in
both the $x$ and $y$ directions by fusing with the $\mathfrak{e}^{z}$
condensate. In the magnetic sector, the $\mathfrak{e}^{z}$ condensate leads to
confinement of fractons. Consider first fracton membrane operators lying in
the $xy$-plane. Such operators do not commute with string operators of the
condensate $\mathfrak{e}^{z}$ lineons, and thus acquire a surface tension.
Naively, $xz$\- and $yz$-plane fracton membrane operators do not acquire any
tension, and indeed they do not feel the $\mathfrak{e}^{z}$ condensate.
However, in the presence of the condensate, the edges of such membranes
running along the $z$-direction become locally detectable, via processes that
create _e.g._ a pair of $\mathfrak{e}^{x}$ lineons, and move one of them in a
small loop encircling the edge before annihilating the lineons to vacuum (the
movement around the loop creates no additional excitations due to the presence
of the $\mathfrak{e}^{z}$ condensate). Nothing prevents adding terms to the
Hamiltonian that effect such processes, giving a line tension to the
$z$-direction edges of the $xz$\- and $yz$-membrane operators. These effects
result in confinement of $\mathfrak{d}^{x,y}_{\mathfrak{m}}$ dipoles, while
leaving $\mathfrak{d}^{z}_{\mathfrak{m}}$ dipoles deconfined. The excitation
content is thus the same as a stack of 2d $\mathbb{Z}_{N}$ gauge theories
oriented along the ${\mathbf{\hat{z}}}$ direction. Indeed, we will see
momentarily that this is the correct identification.

When we further condense a second species of lineon, all of the lineons become
condensed due to the fusion rule
$\mathfrak{e}^{x}\mathfrak{e}^{y}\mathfrak{e}^{z}=1$, and consequently the
magnetic sector is entirely confined. We therefore obtain a trivial
paramagnet.

Figure 12:  The dual representation of the operator $X_{\langle ij\rangle}$ in
terms of a gauge field $\mathcal{Z}_{\langle ij\rangle}$ and two matter spins
$\mathsf{z}_{i},\mathsf{z}_{j}$. The arrows denote the orientations of their
parent links.  Figure 13: An illustration of the operator $B_{p}$ appearing in
$H_{con;z}$ of (67).

Now we will see how these expectations are borne out at the level of explicit
Hamiltonians. To condense lineons we may simply add the kinetic energy term

$H_{\lambda}=-\sum_{a}\lambda^{a}\sum_{l{\,\parallel\,}{\mathbf{\hat{a}}}}(X_{l}+X^{\dagger}_{l}).$
(66)

To identify the phase that results upon condensing $\mathfrak{e}^{z}$ lineons,
we take only $\lambda^{z}\neq 0$ with $\lambda^{x}=\lambda^{y}=0$. Moreover,
we drop the $A^{x}_{i}$ and $A^{y}_{i}$ terms in $H_{XC}$, which tend to gap
out $\mathfrak{e}^{z}$ lineons and do not commute with the $\lambda^{z}$
kinetic energy term. We thus arrive at

$\displaystyle H_{con;z}$ $\displaystyle=$
$\displaystyle-g\sum_{i}A^{z}_{i}-K\sum_{c}B_{c}-\lambda^{z}\sum_{l{\,\parallel\,}{\mathbf{\hat{z}}}}X_{l}$
(67) $\displaystyle-$ $\displaystyle
K^{\prime}\sum_{p{\,\parallel\,}{\mathbf{\hat{z}}}}B_{p}+h.c.$

Here we have also added a term with coefficient $K^{\prime}$, where $B_{p}$ is
a product of $X_{l}$ and $X_{l}^{\dagger}$ over the four links in the
$xy$-plane plaquette $p$, as shown in Fig. 13. This term corresponds precisely
to a process where a pair of $\mathfrak{e}^{x}$ lineons are created, one is
moved around a small $xy$-plane loop, and then the pair is annihilated. As
discussed above, such a process detects the edges of fracton membrane opeators
parallel to ${\mathbf{\hat{z}}}$, and so the $K^{\prime}$ term gives these
edges a line tension. Moreover, if the $K^{\prime}$-term is not added, then an
undesirable unstable ground state degeneracy remains.

To analyze $H_{con;z}$, we first observe that all the terms included commute
with one another. We can then set $X_{l}=1$ for
$l{\,\parallel\,}{\mathbf{\hat{z}}}$, upon which
$B_{c}=B_{p}B_{p^{\prime}}^{\dagger}$, where $p$ (resp. $p^{\prime}$) is the
bottom (resp. top) face of $c$. We see that the $K$-term is redundant and can
be dropped in favor of the $K^{\prime}$-term, giving a stack of
$\mathbb{Z}_{N}$ toric codes as expected.

Now suppose we condense both $\mathfrak{e}^{z}$ and $\mathfrak{e}^{x}$
lineons. This is achieved by the following Hamiltonian:

$H_{con;zx}=-K\sum_{c}B_{c}-\lambda^{z}\sum_{l{\,\parallel\,}{\mathbf{\hat{z}}}}X_{l}-\lambda^{x}\sum_{l{\,\parallel\,}{\mathbf{\hat{x}}}}X_{l}\text{.}$
(68)

Here we have completely dropped the $g$-term to allow both species of lineons
to condense. Setting $X_{l}=1$ for
$l{\,\parallel\,}{\mathbf{\hat{x}}},{\mathbf{\hat{z}}}$ reduces $B_{c}$ to a
product of $X_{l}$ over the four links along ${\mathbf{\hat{y}}}$. The
resulting model is thus a stack of 2d Ising plaquette models (see _e.g._ Ref.
Vijay et al., 2016) along the ${\mathbf{\hat{y}}}$-direction. There is a
ground state degeneracy associated with breaking of subsystem symmetries, but
since we are not imposing these symmetries microscopically, the degeneracy is
not robust. Indeed, we can resolve the degeneracy by adding the term
$-\lambda^{y}\sum_{l{\,\parallel\,}{\mathbf{\hat{y}}}}(X_{l}+h.c)$, which
results in a trivial gapped phase with a unique ground state. Not
coincidentally, this term is the kinetic energy for $\mathfrak{e}^{y}$
lineons, which are automatically condensed, given the fusion rule
$\mathfrak{e}^{x}\mathfrak{e}^{y}\mathfrak{e}^{z}=1$.

To discuss the lineon condensation phase transitions, it is convenient to work
in a dual description in terms of $\mathbb{Z}_{N}$ spins more directly
associated to the lineon excitations, following a logic similar to that used
when discussing condensation of fracton and lineon dipoles. We will begin by
fractionalizing $X_{\langle ij\rangle}$ in terms of dual spins that create
$\mathfrak{e}^{a}$ lineons.121212See Ref. Pretko et al., 2020b for a similar
type of duality. We do this by writing, for each link ${\langle
ij\rangle}{\,\parallel\,}{\mathbf{\hat{a}}}$,

$X_{\langle ij\rangle}=\mathsf{z}^{a}_{i}\mathcal{Z}_{{\langle
ij\rangle}}(\mathsf{z}_{j}^{a})^{\dagger},\qquad Z_{\langle
ij\rangle}=\mathcal{X}_{\langle ij\rangle},$ (69)

where $\mathsf{z}^{a}_{i}$ is to be thought of as creating an
$\mathfrak{e}^{a}$ lineon on site $i$, and where the $\mathcal{Z}_{\langle
ij\rangle},\mathcal{X}_{\langle ij\rangle}$ are dual gauge field variables.
The dual representation of $X_{\langle ij\rangle}$ is shown in figure 12. As
in section IV, in accordance with the lineon fusion rules, at each site we
impose the constraint

$\mathsf{z}_{i}^{x}\mathsf{z}_{i}^{y}\mathsf{z}_{i}^{z}=1.$ (70)

The dual representation and the constraint (70) are invariant under the
$\mathbb{Z}_{N}$ gauge transformations

$\displaystyle\mathsf{z}^{a}_{i}$
$\displaystyle\mapsto\zeta_{i}^{b}(\zeta_{i}^{c})^{*}\mathsf{z}^{a}_{i},$ (71)
$\displaystyle\mathcal{Z}_{\langle
ij\rangle{\,\parallel\,}{\mathbf{\hat{a}}}}$
$\displaystyle\mapsto(\zeta_{i}^{b})^{*}\zeta_{i}^{c}\mathcal{Z}_{\langle
ij\rangle}\zeta_{j}^{b}(\zeta_{j}^{c})^{*},$

where the $\zeta_{i}^{a}$ are valued in the $N$th roots of unity and $a,b,c$
are all distinct. The Gauss laws originating from this gauge redundancy are,
for each $a\neq b$,

$\mathsf{x}^{a}_{i}(\mathsf{x}^{b}_{i})^{\dagger}=\mathcal{X}_{\langle
i,i+{\mathbf{\hat{a}}}\rangle}\mathcal{X}^{\dagger}_{\langle
i,i-{\mathbf{\hat{a}}}\rangle}\mathcal{X}_{\langle
i,i+{\mathbf{\hat{b}}}\rangle}^{\dagger}\mathcal{X}_{\langle
i,i-{\mathbf{\hat{b}}}\rangle},$ (72)

where the terms on the RHS are precisely the duals of the $A^{a}_{i}$
operators appearing in $H_{XC}$ (c.f. figure 2).

There are two types of gauge-invariant operators that create lineons which
will be important in what follows. The first is the dual representation of the
$X_{\langle ij\rangle}$ operators, which hop $\mathfrak{e}^{a}$ lineons along
the ${\mathbf{\hat{a}}}$ direction. The next-simplest term is a “ring-
exchange” operator

$\displaystyle R_{i}^{a}$
$\displaystyle=\mathsf{z}^{a}_{i}\mathcal{Z}_{\langle
i,i+{\mathbf{\hat{b}}}\rangle}(\mathsf{z}^{a}_{i+{\mathbf{\hat{b}}}})^{\dagger}\mathcal{Z}_{\langle
i+{\mathbf{\hat{b}}},i+{\mathbf{\hat{b}}}+{\mathbf{\hat{c}}}\rangle}^{\dagger}\mathsf{z}^{a}_{i+{\mathbf{\hat{b}}}+{\mathbf{\hat{c}}}}$
(73) $\displaystyle\qquad\times\mathcal{Z}^{\dagger}_{\langle
i+{\mathbf{\hat{c}}},i+{\mathbf{\hat{b}}}+{\mathbf{\hat{c}}}\rangle}(\mathsf{z}^{a}_{i+{\mathbf{\hat{c}}}})^{\dagger}\mathcal{Z}_{\langle
i,i+{\mathbf{\hat{c}}}\rangle},$

where as usual $a,b,c$ are all distinct. This operator creates four
$\mathfrak{e}^{a}$ lineons at the corners of a plaquette normal to
${\mathbf{\hat{a}}}$ and with a corner at the site $i$, and will be seen to
play an important role in our analysis of the condensation transition.

The Hamiltonian can now be written as

$\displaystyle H_{con}$ $\displaystyle=-g\sum_{a\neq
b}(\mathsf{x}_{i}^{a})^{\dagger}\mathsf{x}_{i}^{b}-K\sum_{c}\mathcal{B}_{c}$
(74) $\displaystyle\qquad-\sum_{a}\lambda^{a}\sum_{{\langle
ij\rangle}{\,\parallel\,}{\mathbf{\hat{a}}}}\mathsf{z}^{a}_{i}\mathcal{Z}_{\langle
ij\rangle}(\mathsf{z}^{a}_{j})^{\dagger}+h.c,$

where $\mathcal{B}_{c}$ is $B_{c}$ but with $\mathcal{Z}_{\langle ij\rangle}$
operators instead of $X_{\langle ij\rangle}$ operators. Magnetic excitations
are gapped, so as in the previous examples we can ignore fluctuations of the
$\mathcal{Z}_{\langle ij\rangle}$ generalized gauge field, which indeed are
not present in (74).

If we set $\mathcal{Z}_{\langle ij\rangle}=1$, we obtain an effective matter
theory for the lineons, which is valid when considering correlation functions
of local, gauge-invariant operators. This effective matter theory enjoys
planar $\mathbb{Z}_{N}$ subsystem symmetries, arising from the planar
$\mathbb{Z}_{N}$ conservation laws for lineons in the X-cube phase discussed
in Sec. II. Concretely, the Hamiltonian and any perturbations are required to
be invariant under the transformation

$\mathsf{z}^{a}_{i}\to\zeta^{b}_{i}\zeta^{c}_{i}\mathsf{z}^{a}_{i}\text{,}$
(75)

where $a,b,c$ are all distinct, and where $\zeta^{b}_{i}$ and $\zeta^{c}_{i}$
are constant within planes normal to ${\mathbf{\hat{b}}}$ and
${\mathbf{\hat{c}}}$, respectively.

###  V.2 Phase transitions for $N\leq 4$

We now examine the nature of the phase transitions that occur when the lineons
condense. We will first consider only the simplest case when a single species
of lineon condenses, which as in the previous section we will take to be
$\mathfrak{e}^{z}$. We will therefore take $\lambda^{a}=(0,0,\lambda)$, and
consider the condensation transition that occurs as a function of $g/\lambda$.

In this limit, $\mathsf{x}^{x},\mathsf{x}^{y}$ both commute with the
Hamiltonian (74), and to study the phase transition we may work in the
$\mathsf{x}^{x}=\mathsf{x}^{y}=1$ subspace. In this subspace the Hamiltonian
(74) clearly reduces to an array of 1d $\mathbb{Z}_{N}$ clock model wires
oriented along the ${\mathbf{\hat{z}}}$ direction, all of which are coupled to
the (non-fluctuating) $\mathbb{Z}_{N}$ gauge field $\mathcal{Z}_{\langle
ij\rangle}$. In this limit, the critical point is simply described by a
decoupled array of 1d critical clock models, with the critical fluctuations
occurring on each wire independently.

Of course the limit where the $\mathbb{Z}_{N}$ chains are completely decoupled
is very fine-tuned, and we are interested in the fate of the decoupled fixed
point under generic perturbations respecting the effective subsystem symmetry.
This can be determined by performing an RG analysis – if all of the
perturbations which couple different chains are irrelevant, the decoupled
fixed point will correctly describe the phase transition in a finite region of
parameter space. If some of the perturbations are relevant, they may drive the
system into a nontrivial coupled fixed point, or may lead to the transition
being driven first order. Using our knowledge of the critical properties of
the 1d $\mathbb{Z}_{N}$ clock model, the RG analysis can be performed
straightforwardly. Two classes of terms need to be considered, those which are
a) invariant under separate $\mathbb{Z}_{N}$ transformations on each wire or
are b) built from the ring exchange terms of (73).

####  V.2.1 $N=2$

First, consider the case where $N=2$. The decoupled critical point is then
described by an array of critical Ising CFTs. Energy-energy interactions are
an important class of couplings between pairs of chains. Labeling the chains
by $w$ and assuming translation symmetry, these couplings give us a
perturbation to the decoupled fixed point which can be written as

$\delta S_{\varepsilon}=\sum_{w\neq w^{\prime}}C_{w-w^{\prime}}\int
dz\,d\tau\,\varepsilon_{w}\varepsilon_{w^{\prime}},$ (76)

where $\varepsilon_{w}$ is the energy operator on the chain $w$ (note that
since there is no $\varepsilon^{2}$ operator in the 1d Ising CFT, we only sum
over $w^{\prime}\neq w$). Since the critical fluctuations are only in the
$z$-$\tau$ spacetime plane, the relevance of $\delta S$ is obtained by
comparing the scaling dimension of $\varepsilon_{w}\varepsilon_{w^{\prime}}$
to 2. Since the dimension of the energy operator in the 1d Ising model is 1,
the couplings $C_{w,w^{\prime}}$ are all marginal at tree level.

However, it is unnecessary to actually compute the beta functions for the
$C_{w-w^{\prime}}$ couplings. This is because there is another gauge-invariant
coupling which is relevant, which as we will see destabilizes all of the
decoupled critical points when $N\leq 4$. This is the ring-exchange term (73),
which at the critical point is written as (dropping the $\mathcal{Z}$ gauge
fields as they do not affect calculations of scaling dimensions at the
critical point)

$\displaystyle\delta S_{r}$ $\displaystyle=r\sum_{i,a}\int\int
dz\,d\tau\,R_{i}^{a}$ (77) $\displaystyle\rightarrow
r\sum_{w_{1,\dots,4}\in\square}\int
dz\,d\tau\,\sigma_{w_{1}}\sigma_{w_{2}}\sigma_{w_{3}}\sigma_{w_{4}},$

where $\sigma_{w}$ is the spin operator on wire $w$ and the notation
$w_{1,\dots,4}\in\square$ means a sum over configurations of four wires
$w_{1,\dots,4}$ whose $x,y$ coordinates form the corners of a square with unit
normal along ${\mathbf{\hat{z}}}$. Since the scaling dimension of $\sigma$ is
$\Delta_{\sigma}=1/8$, the scaling dimension of the above term is

$\Delta_{r}=4\Delta_{\sigma}=\frac{1}{2}<2,$ (78)

which is very relevant. Higher-order contributions to the beta function for
$r$ do not lead to fixed points accessible at small couplings. We therefore
conclude that the transition is likely generically first-order, in agreement
with suggestions from numerics.Slagle and Kim (2017a) The same conclusion
holds when mutiple species of lineons condense simultaneously.

####  V.2.2 $N=3$

Now consider $N=3$. In this case the putative decoupled critical point
corresponds to an array of $\mathbb{Z}_{3}$ clock models. At the critical
point of the $\mathbb{Z}_{3}$ chain, the scaling dimension of the spin
operator $\sigma$ is $\Delta_{\sigma}=2/15$, while that of the energy operator
is $\Delta_{\varepsilon}=4/5$.Dotsenko (1984) Therefore both the ring-exchange
term and the energy-energy couplings are relevant, and the transitions are
again expected to be first order.

####  V.2.3 $N=4$

Next up is $N=4$. While in this case there is a one-parameter critical line of
fixed points (the Ashkin-Teller line), the scaling dimension of the spin
operator on each chain is always $\Delta_{\sigma}=1/8$.Delfino and Grinza
(2004) As such the ring-exchange term is again relevant, and destabilizes the
decoupled critical point.

###  V.3 Phase transitions for $N>4$

For $N>4$, the 1d chains are more conveniently dealt with using a continuum XY
description.José et al. (1977) We do this by writing the $\mathsf{z}^{z}$ and
$\mathsf{x}^{z}$ variables as exponentials of slowly-varying fields
$\Phi_{w}(z)$ and $\Theta_{w}(z)$, so that in the IR near the phase transition
we have the approximate identifications

$\displaystyle\mathsf{z}^{z}_{(w,z)}$
$\displaystyle\sim\exp\left(i\Phi_{w}(z)\right),$ (79)
$\displaystyle\mathsf{x}^{z}_{(w,z)}$
$\displaystyle\sim\exp\left(i\frac{1}{N}\partial_{z}\Theta_{w}(z)\right)=\exp\left(i\frac{{2\pi}}{N}\pi_{\Phi_{w}}(z)\right),$

where $\pi_{\Phi_{w}}=\frac{1}{{2\pi}}\partial_{z}\Theta_{w}$ is the momentum
conjugate to $\Phi_{w}$, and where a suitable continuum limit is taken in the
${\mathbf{\hat{z}}}$ direction, assuming a continuous transition. The
$\mathbb{Z}_{N}$ nature of the problem can be accounted for by including a
$\mathbb{Z}_{N}$ anisotropy term $\cos(N\Phi_{w})$ in the action.

The most general action for the putative fixed point can then be written as
$S=S_{0}+S_{I}$, where $S_{I}$ contains interactions and where the first term
represents a non-interacting quadratic fixed point:

$\displaystyle S_{0}$ $\displaystyle=\frac{1}{{4\pi}}\sum_{w,w^{\prime}}\int
dz\,d\tau\,\big{(}2\delta_{w,w^{\prime}}\partial_{\tau}\Phi_{w}\partial_{z}\Theta_{w^{\prime}}$
(80)
$\displaystyle\qquad+g_{w-w^{\prime}}\partial_{z}\Theta_{w}\partial_{z}\Theta_{w^{\prime}}+\eta_{w-w^{\prime}}\partial_{z}\Phi_{w}\partial_{z}\Phi_{w^{\prime}}\big{)},$

where we have assumed translation invariance for the interchain derivative
couplings $g_{w-w^{\prime}},\eta_{w-w^{\prime}}$.

$S_{I}$ contains the $\mathbb{Z}_{N}$ anisotropy and further cosines in
$\Phi_{w},\Theta_{w}$ allowed by gauge invariance, with the most important
ones for the present problem being

$\displaystyle S_{I}$ $\displaystyle=a^{-2}\int
dz\,d\tau\,\Big{[}u\sum_{w}\cos(N\Phi_{w})$ (81)
$\displaystyle\qquad+r\sum_{w_{1},\dots,w_{4}\in\square}\cos\Big{(}\Phi_{w_{1}}-\Phi_{w_{2}}+\Phi_{w_{3}}-\Phi_{w_{4}}\Big{)}$
$\displaystyle\qquad+s\sum_{w}\cos(\Theta_{w})+\cdots\Big{]},$

where $a$ is the lattice spacing. Despite the fact that we naively need
$u\rightarrow\infty$ to enforce the $\mathbb{Z}_{N}$ nature of the problem,
$u$ can in fact be treated perturbatively.131313This can be shown by requiring
that the known physics of the XY model be recovered in the
$N\rightarrow\infty$ limit; see Refs. Elitzur et al., 1979; Fröhlich and
Spencer, 1981.

####  V.3.1 Decoupled limit

Let us first compute the scaling dimensions of the terms in $S_{I}$ in the
case where there are no interchain derivative couplings, i.e. let us first
consider setting

$g_{w-w^{\prime}}=g\delta_{w,w^{\prime}},\qquad\eta_{w-w^{\prime}}=\eta\delta_{w,w^{\prime}}.$
(82)

The scaling dimensions of the most relevant gauge-invariant cosines are then
found to be

$\displaystyle\Delta_{u}$ $\displaystyle=\frac{N^{2}}{2R^{2}}$ (83)
$\displaystyle\Delta_{r}$ $\displaystyle=\frac{2}{R^{2}}$
$\displaystyle\Delta_{s}$ $\displaystyle=\frac{R^{2}}{2},$

where we have defined

$R^{2}\equiv\sqrt{\eta/g}.$ (84)

The relevance of a given term in $S_{I}$ is determined by comparing the
appropriate scaling dimension to two. When the field term is small, so that
the lineons are not condensed, $R^{2}$ is also small. Here the
$\cos(\Theta_{w})$ operators are relevant, pinning the values of the
$\Theta_{w}$ fields. This regime occurs for $R^{2}<4$. In the limit of large
field where the lineons are condensed, $R^{2}$ is large, and here the
$\cos(N\Phi_{w})$ term is relevant. This happens for $R^{2}>N^{2}/4$. In the
absence of the ring exchange term this would give an intermediate massless
regime for $4<R^{2}<N^{2}/4$ with no relevant perturbations. However, we see
that the ring exchange term is relevant for all $R^{2}>1$, meaning that there
is in fact no choice of $R^{2}$ for which all three terms in (83) are
irrelevant.

It remains to understand the nature of the intermediate $R^{2}$ regime when
the ring exchange term is the dominant relevant perturbation. What happens in
this case is explained in appendix C. In brief, we find that the ring-exchange
term leads to a putative fixed point corresponding to an anisotropic version
of the exciton Bose liquid,Paramekanti et al. (2002); You et al. (2020a);
Seiberg and Shao (2020a, b) but that the three-dimensional nature of the
problem produces an instability with respect to lineon dipole condensation. We
then conclude that when the ring exchange term dominates, the transition is
expected be first order.

####  V.3.2 General derivative couplings

That the decoupled fixed point is destabilized by the ring-exchange terms does
not necessarily mean that there do not exist stable fixed points described by
the quadratic action (80) — it only means that any such fixed point must owe
its existence to interchain derivative couplings which frustrate ring-exchange
processes. In this section we show numerically that a relatively simple choice
for the derivative couplings can be made which renders the Gaussian fixed
point stable with respect to low-body finite-range cosines. Whether or not the
chosen couplings ensure stability with respect to all possible deformations is
a more complicated question, which we leave to future work. However, we also
argue that there exist other, more complicated, choices of interchain
couplings that result in a truly stable fixed point.

We begin by Fourier transforming in the $x$,$y$ directions normal to the
chains, writing $S_{0}$ as

$\displaystyle S_{0}$
$\displaystyle=\frac{1}{{4\pi}}\int\frac{d^{2}\mathbf{q}}{({2\pi})^{2}}\,\int
dz\,d\tau\,\Bigg{[}2\partial_{\tau}\Phi_{\mathbf{q}}^{*}\partial_{z}\Theta_{\mathbf{q}}$
(85) $\displaystyle\qquad
v_{\mathbf{q}}\left(\frac{1}{R^{2}_{\mathbf{q}}}\partial_{z}\Theta_{\mathbf{q}}^{*}\partial_{z}\Theta_{\mathbf{q}}+R^{2}_{\mathbf{q}}\partial_{z}\Phi_{\mathbf{q}}^{*}\partial_{z}\Phi_{\mathbf{q}}\right)\Bigg{]},$

where the $q_{x}$ and $q_{y}$ momentum integrations run over $[0,2\pi)$ (with
the interchain spacing set to unity), and where we have defined

$v_{\mathbf{q}}\equiv\sqrt{\eta_{\mathbf{q}}g_{\mathbf{q}}},\qquad
R^{2}_{\mathbf{q}}\equiv\sqrt{\eta_{\mathbf{q}}/g_{\mathbf{q}}}.$ (86)

For stability we require that the couplings be chosen so that
$0<R_{\mathbf{q}}^{2}<\infty$ for all $\mathbf{q}$; in what follows we will
assume that this is the case.

The scaling dimensions of perturbations to this fixed point can be computed as
various integrals of $R^{2}_{\mathbf{q}}$-dependent functions. Therefore to
demonstrate the existence of a stable fixed point, all that remains is to find
an appropriate choice of $R^{2}_{\mathbf{q}}$ which renders all the allowed
perturbations irrelevant. The problem of finding such a function
$R^{2}_{\mathbf{q}}$ is very closely related to the problem of constructing
examples of “sliding Luttinger liquids”,Mukhopadhyay et al. (2001); Vishwanath
and Carpentier (2001) which are anisotropic metallic phases built from
strongly coupled arrays of Luttinger liquids.

At the fixed point (85), general cosines of the $\Phi$ and $\Theta$ variables
have scaling dimensions

$\displaystyle\Delta_{\cos(\sum_{w}\alpha_{w}\Phi_{w})}$
$\displaystyle=\frac{1}{2}\int\frac{d^{2}\mathbf{q}}{({2\pi})^{2}}\frac{|\alpha_{\mathbf{q}}|^{2}}{R_{\mathbf{q}}^{2}}$
(87) $\displaystyle\Delta_{\cos(\sum_{w}\beta_{w}\Theta_{w})}$
$\displaystyle=\frac{1}{2}\int\frac{d^{2}\mathbf{q}}{({2\pi})^{2}}R_{\mathbf{q}}^{2}|\beta_{\mathbf{q}}|^{2},$

where $\alpha_{w},\beta_{w}$ are integer-valued functions of the wire index
$w$, which are non-zero only for finitely many $w$.

Figure 14:  An illustration of the geometry of a ring exchange term
$R^{\mathbf{a}\mathbf{b}}$. Here the vectors $\mathbf{a},\mathbf{b}$ are in
the $x$-$y$ plane and connect four wires to form a rectangle as shown.

The simplest gauge-invariant cosines are $\cos(N\Phi_{w})$,
$\cos(\Theta_{w})$, $\cos(\Theta_{w}-\Theta_{w+\mathbf{a}})$, and
$\cos(R^{ab})$, where now the more general ring-exchange operator
$R^{\mathbf{a}\mathbf{b}}$ is defined as (again omitting the unimportant
$\mathcal{Z}$ gauge field)

$R^{\mathbf{a}\mathbf{b}}=\Phi_{w}-\Phi_{w+\mathbf{a}}+\Phi_{w+\mathbf{a}+\mathbf{b}}-\Phi_{w+\mathbf{b}},$
(88)

where $\mathbf{a}\perp\mathbf{b}$ form the edges of a rectangle (see figure
14). Using the above formula, we see that these have the scaling dimensions

$\displaystyle\Delta_{\cos(N\Phi_{w})}$
$\displaystyle=\frac{N^{2}}{2R^{2}_{\mathbf{0}}}$ (89)
$\displaystyle\Delta_{\cos(\Theta_{w})}$
$\displaystyle=\frac{R^{2}_{\mathbf{0}}}{2}$
$\displaystyle\Delta_{\cos(\Theta_{w}-\Theta_{w+\mathbf{a}})}$
$\displaystyle=\frac{1}{2}\int\frac{d^{2}\mathbf{q}}{({2\pi})^{2}}R^{2}_{\mathbf{q}}(1-\cos(\mathbf{a}\cdot\mathbf{q}))$
$\displaystyle\Delta_{\cos(R^{\mathbf{a}\mathbf{b}})}$
$\displaystyle=\int\frac{d^{2}\mathbf{q}}{({2\pi})^{2}}\frac{|1-e^{-i\mathbf{a}\cdot\mathbf{q}}-e^{-i\mathbf{b}\cdot\mathbf{q}}+e^{i(\mathbf{a}+\mathbf{b})\cdot\mathbf{q}}|^{2}}{2R^{2}_{\mathbf{q}}}.$

By taking $N$ large enough we see that we can always make the first term (the
$\mathbb{Z}_{N}$ anisotropy) irrelevant; hence we can ignore it for the time
being. For the remaining terms, we see that in order to make the above scaling
dimensions large we should take $R_{\mathbf{q}}^{2}$ to be a) everywhere
nonzero, b) small at the wavevectors where the numerator in the integrand of
the expression for $\Delta_{\cos(R^{\mathbf{a}\mathbf{b}}_{\Phi})}$ is small,
and c) large enough at other wavevectors so that the cosines of $\Theta$ are
kept irrelevant.

We will not concern ourselves with a completely general analysis of which
functional forms of $R^{2}_{\mathbf{q}}$ do the job, and will be content with
simply demonstrating that one particular choice of $R^{2}_{\mathbf{q}}$ works.

One such choice, which was employed in Ref. Mukhopadhyay et al., 2001 in the
context of a slightly different problem, is

$R^{2}_{\mathbf{q}}=\eta_{0}\left(1+\eta_{1}[\cos(q_{x})+\cos(q_{y})]+\eta_{2}\cos(q_{x})\cos(q_{y})\right)^{2},$
(90)

with $\eta_{0}>0$. For stability $R_{\mathbf{q}}^{2}$ must be nonzero for all
$\mathbf{q}$, which ends up restricting $\eta_{1},\eta_{2}$ to be such that

$-\frac{1}{2}<\eta_{1}<\frac{1}{2},\qquad 2|\eta_{1}|-1<\eta_{2}<1.$ (91)

The dimensions of the first two operators in (89) are easy to compute, and are

$\displaystyle\Delta_{\cos(\Theta)}$
$\displaystyle=\frac{\eta_{0}}{2}\left(1+\eta_{1}^{2}+\frac{\eta_{2}^{2}}{4}\right),$
(92) $\displaystyle\Delta_{\cos(N\Phi)}$
$\displaystyle=\frac{N^{2}}{2\eta_{0}(1+\eta_{1}^{2}+\eta_{2}^{2}/4)}.$

Figure 15:  A region of stability (shaded gray area) in the
$\eta_{1}$-$\eta_{2}$ plane, where the fixed-point model has no relevant
perturbations. In the upper region marked “$\exists$ rel couplings” there
exist relevant deformations to the fixed point, while in the lower region the
model is unstable, as there is always some $\mathbf{q}$ for which
$R^{2}_{\mathbf{q}}=0$.

To figure out if there are any regions in the $\eta_{1}$-$\eta_{2}$ plane
where all of the scaling dimensions in (89) are greater than 2, we resort to a
numerical search. For these purposes we restrict the search to couplings
involving fields up to 10th nearest neighbors. We find small regions of
stability near the boundaries of the region in which $R^{2}_{\mathbf{q}}$ is
strictly positive; one such region is shown as the gray shaded area in figure
15. In this region $\cos(N\Phi)$ is irrelevant for all $N>4$. Note that
although the fixed points described here are very anisotropic, they are
nevertheless not completely decoupled (indeed, the coupling is crucial for
their stability), unlike the fixed points considered in section III.3.

Of course, demonstrating that the scaling dimensions in (89) are all greater
than 2 does not actually guarantee stability, as there may be more complicated
higher-body / longer-range operators which are relevant. Note that as a matter
of principle it is always possible to construct a stable fixed point with
sufficiently complicated interchain interactions, as was shown in Ref.
Plamadeala et al., 2014. In that work, authors showed that a single quantum
wire of interacting bosons has a stable gapless “perfect metal” phase, where
all cosine operators are irrelevant, when the number of channels is $8k$, with
$k\geq 3$. The action is quadratic and crucially involves interchannel
derivative couplings, as in $S_{0}$ of (80). It was also pointed out that
arrays of such perfect metals give rise to stable gapless phases in 2d and 3d.
In the present case, we can employ these results by bundling the wires of our
model into groups of 24 (which breaks full translation invariance), and
choosing interchain couplings within each bundle to put it in a perfect metal
phase. The resulting system of decoupled bundles is then a stable gapless
phase.

Even with the present choice of relatively simple interchain couplings, any
relevant cosines must necessarily involve rather high-body
interactions141414The scaling dimensions of higher-body interactions of the
$\Theta$ fields are easy to check, since $R_{\mathbf{q}}^{2}$ only contains a
handful of Fourier modes. $1/R_{\mathbf{q}}^{2}$ on the other hand contains an
infinite number of Fourier modes, which makes a systematic search of the
scaling dimensions of the $\Phi$ cosines more difficult. and / or involve very
long-range interchain couplings. Therefore as in the case of the “almost
perfect metals” of Ref. Murthy and Nayak, 2020, any relevant deformations are
likely to have very small bare values, meaning that the examples constructed
above will only show instabilities at very low temperatures.

With the aforementioned caveats about stability in mind, we thus arrive at a
scenario where the XC model is separated from the condensed phase by a
massless critical phase extending over a finite region of parameter space. The
transition from the XC phase to the intermediate massless region is likely to
be of BKT type, since the operator which becomes relevant in the XC phase is
the cosine $\cos(\Theta)$. The nature of the transition from the massless
region to the condensed phase is not completely clear, and may depend on the
type of operator which becomes relevant at the transition.

Finally, we note that although in this section we have only considered the
transitions that occur as a single species of lineon is condensed, the
generalization to the case where multiple species condense simultaneously is
straightforward, due to the fact that couplings between perpendicular wires
are always less relevant than those between parallel wires. In particular,
when multiple lineons condense we continue to find no stable fixed points for
$N\leq 4$.

##  VI Discussion

In this paper we have discussed various types of condensation transitions in
the $\mathbb{Z}_{N}$ XC model, with the phases on both sides of the
transitions described in terms of various types of gauge theories. We have
identified continuous condensation phase transitions, and intermediate gapless
phases, provided that $N>4$, although our analysis does not rule out the
existence of continuous transitions for smaller $N$. Clearly there is a lot
more work along these lines that may be done.

One obvious extension of the present work is to consider similar decoupled
critical points in other fracton models, and to understand them using our
existing knowledge of 1d and 2d critical phenomena. This strategy is likely to
yield many other interesting examples.

One interesting question to ask is then whether or not there exists a
continuous phase transition between a fracton phase and a phase without
fracton order, where the critical modes at the transition fluctuate in the
full four-dimensional spacetime. In fact, one such critical point was proposed
in Ref. Vijay, 2017. This work discussed the transition between the XC model
and decoupled stacks of two-dimensional $\mathbb{Z}_{N}$ gauge theories. Based
on a duality between a model of coupled $\mathbb{Z}_{N}$ gauge theory layers
and a 3d $\mathbb{Z}_{N}$ gauge theory, it was proposed that this transition
is first order for $N\leq 4$, but that there may be a continuous transition
out of the XC phase for $N>4$. This transition is dual to that between the
deconfined phase of 3d $\mathbb{Z}_{N}$ gauge theory and a massless Coulomb
phase, which was claimed to be continuous in early Monte Carlo studies.Creutz
et al. (1979) However, we believe that this transition is in fact likely
weakly first order. Indeed, starting from the Coulomb phase, one may obtain
the $\mathbb{Z}_{N}$ gauge theory by condensing electric particles with charge
$N$. If we assume a second order transition, the critical point can presumably
be described by a charge-$N$ Higgs field coupled to a $U(1)$ gauge field. The
Higgs transition with $n\lesssim 365$ flavors in 3d is however generically
made first order by fluctuations,Halperin et al. (1974) and hence the
transition from the XC model to the intermediate massless phase is also likely
to be weakly first order.

One possible route to a continuous transition with 3d character could
potentially lie in Higgs transitions into the XC phase, of the type studied in
Refs.Ma et al. (2018); Bulmash and Barkeshli (2018) The nonstandard dispersion
of the matter fields in these examples may help to stabilize against a
fluctuation-induced first-order transition, although the ultimate character of
the phase transition may also end up being quasi 2d. Another possible route
lies in finding a theory where the decoupled fixed point is unstable, but can
be shown to flow to a stable fixed point with nontrivial couplings between
different planes / chains. We leave a more detailed treatment of these
possibilities to future work.

## Acknowledgements

EL thanks Dave Aasen, Han Ma, T. Senthil, and Yizhi You for discussions. EL is
supported by the Fannie and John Hertz Foundation and the NDSEG fellowship.
The research of MH is supported by the U.S. Department of Energy, Office of
Science, Basic Energy Sciences (BES) under Award number DE-SC0014415. This
work was also partly supported by the Simons Collaboration on Ultra-Quantum
Matter, which is a grant from the Simons Foundation (651440, MH).

##  Appendix A $\mathbb{Z}_{N}$ fracton dipole condensation

Figure 16:  A reference for the pattern of Hermitian conjugations appearing in
the Gauss law (93). Links shaded in orange denote $X_{\mathsf{l}}^{\dagger}$
operators, and those shaded in teal denote $X_{\mathsf{l}}$ operators. Each
orange plaquette denotes a $\mathcal{X}^{\dagger}_{\mathsf{p}}$ operator, and
the unshaded plaquettes denote $\mathcal{X}_{\mathsf{p}}$ operators.  Figure
17:  The kinetic terms employed when condensing dipoles (plus their Hermitian
conjugates, not shown). The standard labeling of $l_{1,\dots,4}$ for each
orientation of plaquette is indicated in the left column.

In this appendix we discuss how the generalized gauge theory employed in the
discussion of fracton dipole condensation in section III.1 is generalized to
the $\mathbb{Z}_{N}$ case. The generalization is rather straightforward, with
the difficulties lying only in keeping track of the correct ways to take
Hermitian conjugates in various expressions.

As in section III.1 we dualize the XC model by writing it in terms of
$\mathbb{Z}_{N}$ matter qubits $\mathsf{x},\mathsf{z}$ on the vertices
$\mathsf{i}$ of the dual lattice, and $\mathbb{Z}_{N}$ gauge qubits
$\mathcal{X},\mathcal{Z}$ on the plaquettes $\mathsf{p}$ of the dual lattice.

We then add additional $\mathbb{Z}_{N}$ qubits $X,Z$ on the links of dual
lattice which represent the dipole degrees of freedom. The gauss law reads, in
analogy to (8),

$\mathsf{x}_{\mathsf{i}}=\prod_{\mathsf{l},\mathsf{p}\in\partial\mathsf{i}_{out}}\mathcal{X}_{\mathsf{p}}^{\dagger}X_{\mathsf{l}}^{\dagger}\prod_{\mathsf{l}^{\prime},\mathsf{p}^{\prime}\in\partial\mathsf{i}_{in}}\mathcal{X}_{\mathsf{p}^{\prime}}X_{\mathsf{l}^{\prime}}.$
(93)

Here $\partial\mathsf{i}_{out}$ consists of the collection of outward-oriented
plaquettes and links neihboring $\mathsf{i}$, while $\partial\mathsf{i}_{in}$
consists of the inward-oriented plaquettes and links. Here outward-oriented
(inward-oriented) plaquettes are those whose centers have coordinates with
respect to $\mathsf{i}$ which are all of the same sign (of different signs),
and are indicated as shaded (not shaded) in figure 16. Outward-oriented links
are those oriented parallel to the coordinate axes (orange in figure 16), and
inward-oriented links are those oriented anti-parallel (cyan in figure 16).

With this notation and the above Gauss law, the gauge-invariant kinetic terms
which hop dipoles across a given plaquette are then
$\mathcal{Z}_{\mathsf{p}}Z_{\mathsf{l}_{1}}Z_{\mathsf{l}_{2}}^{\dagger}$ and
$\mathcal{Z}_{\mathsf{p}}Z_{\mathsf{l}_{3}}Z_{\mathsf{l}_{4}}^{\dagger}$,
together with their Hermitian conjugates. The correct $\mathbb{Z}_{N}$
generalization of the Hamiltonian (9) is then

$\displaystyle H_{con}$
$\displaystyle=-g\sum_{\mathsf{i}}\mathsf{x}_{\mathsf{i}}-h\sum_{\mathsf{l}}X_{\mathsf{l}}-\lambda\sum_{\mathsf{p}}\mathcal{Z}_{\mathsf{p}}(Z_{\mathsf{l}_{1}}Z_{\mathsf{l}_{2}}^{\dagger}+Z_{\mathsf{l}_{3}}Z_{\mathsf{l}_{4}}^{\dagger})+h.c,$
(94)

where as before we have ignored the dual representation of the
$A^{a}_{\mathsf{i}}$ terms, which will not be important for discussing the
condensed phase.

To demonstrate that $H_{con}$ above is equivalent to the Hamiltonian for
deconfined $\mathbb{Z}_{N}$ gauge theory in the $h/\lambda\rightarrow 0$
limit, we again need to perform a unitary transformation which decouples the
plaquette degrees of freedom. This is done by directly generalizing the
analysis of the $\mathbb{Z}_{2}$ case.

To begin, define the operators

$\Pi_{\mathsf{p}k}=\frac{1}{N}\sum_{m\in\mathbb{Z}_{N}}(Z_{\mathsf{l}_{1}}Z_{\mathsf{l}_{2}}^{\dagger})^{m}e^{ikm},$
(95)

which form a complete set of Hermitian projectors. Now define the unitaries

$U_{\mathsf{p}}=\sum_{k}\Pi_{\mathsf{p}k}\mathcal{X}^{k}_{\mathsf{p}},$ (96)

and as before let $\mathcal{U}\equiv\prod_{\mathsf{p}}U_{\mathsf{p}}$. The
unitary $\mathcal{U}$ conjugates operators as

$\displaystyle\mathcal{U}^{\dagger}\mathcal{X}_{\mathsf{p}}\mathcal{U}$
$\displaystyle=\mathcal{X}_{\mathsf{p}}$ (97)
$\displaystyle\mathcal{U}^{\dagger}\mathcal{Z}_{\mathsf{p}}\mathcal{U}$
$\displaystyle=Z_{\mathsf{l}_{1}}^{\dagger}Z_{\mathsf{l}_{2}}\mathcal{Z}_{\mathsf{p}}$
$\displaystyle\mathcal{U}^{\dagger}X_{\mathsf{l}_{1}}\mathcal{U}$
$\displaystyle=\mathcal{X}_{\mathsf{p}}X_{\mathsf{l}_{1}}$
$\displaystyle\mathcal{U}^{\dagger}X_{\mathsf{l}_{2}}\mathcal{U}$
$\displaystyle=\mathcal{X}_{\mathsf{p}}^{\dagger}X_{\mathsf{l}_{2}}$
$\displaystyle\mathcal{U}^{\dagger}X_{\mathsf{l}_{3,4}}\mathcal{U}$
$\displaystyle=X_{\mathsf{l}_{3,4}}$
$\displaystyle\mathcal{U}^{\dagger}Z_{\mathsf{l}_{\alpha}}\mathcal{U}$
$\displaystyle=Z_{\mathsf{l}_{\alpha}},$

which can be verified using the orthogonality and completeness of the
$\Pi_{\mathsf{p}k}$ as well as the relations

$\displaystyle X_{\mathsf{l}_{1}}\Pi_{\mathsf{p}k}$
$\displaystyle=\Pi_{\mathsf{p}k-1}X_{\mathsf{l}_{1}}$ (98) $\displaystyle
X_{\mathsf{l}_{2}}\Pi_{\mathsf{p}k}$
$\displaystyle=\Pi_{\mathsf{p}k+1}X_{\mathsf{l}_{2}}$ $\displaystyle
X_{\mathsf{l}_{3,4}}\Pi_{\mathsf{p}k}$
$\displaystyle=\Pi_{\mathsf{p}k}X_{\mathsf{l}_{3,4}}.$

Figure 18:  The plaquette terms appearing in the $\mathbb{Z}_{N}$ gauge theory
(Hermitian conjugates not shown).

Using these relations and referring to the pattern of Hermitian conjugation in
figure 17, one sees that under conjugation by $\mathcal{U}$, the Gauss law
constraint maps to

$\mathsf{x}_{\mathsf{i}}=\prod_{\mathsf{l}\in\partial\mathsf{i}_{out},\mathsf{l}^{\prime}\in\partial\mathsf{i}_{in}}X_{\mathsf{l}}X_{\mathsf{l}^{\prime}}^{\dagger},$
(99)

which is the desired Gauss law of $\mathbb{Z}_{N}$ gauge theory.

By conjugating the kinetic term in $H$ with $\mathcal{U}$, one sees that as in
the $\mathbb{Z}_{2}$ case, the conjugation generates a set of plaquette terms
identical to that of the standard presentation of $\mathbb{Z}_{N}$ gauge
theory, together with a trivial paramagnet coming from the plaquette
variables. After getting rid of the plaquette degrees of freedom, the
Hamiltonian is

$\displaystyle H_{con}^{\prime}$
$\displaystyle=\mathcal{U}^{\dagger}H_{con}\mathcal{U}$ (100)
$\displaystyle=-g\sum_{\mathsf{i}}\mathsf{x}_{\mathsf{i}}-\lambda\sum_{\mathsf{p}}B_{\mathsf{p}}+h.c,$

where the $B_{\mathsf{p}}$ terms are given by the operators appearing in
figure 18. This completes our identification of the condensed phase with
deconfined $\mathbb{Z}_{N}$ gauge theory.

##  Appendix B RG flow for stacked Ising models

In this appendix we argue that $L$ stacked layers of critical 2d Ising models
are always unstable with respect to interlayer interactions that preserve the
$\mathbb{Z}_{2}^{L}$ subsystem symmetry

$\mathbb{Z}_{2}^{L}:\psi_{i}\mapsto f_{i}\psi_{i},\qquad f_{i}\in\\{1,-1\\},$
(101)

where $\psi_{i}$ is the Ising field on layer $i$.

We will address this question within the context of the $\varepsilon$
expansion. Our starting point is then the critical action

$\displaystyle S$ $\displaystyle=\int
d^{3}x\,\Big{(}\,\frac{1}{2}\sum_{i=1}^{L}(\partial\psi_{i})^{2}+\frac{\Lambda}{8}\sum_{i,j=1}^{L}\psi_{i}^{2}\psi_{j}^{2}g_{i,j}\Big{)},$
(102)

where the couplings $g_{i,j}$ are symmetric and dimensionless, $\Lambda$ is a
UV cutoff, and where the absence of mass terms denotes tuning to the critical
point.

Enumerating all the symmetric RG fixed points for even moderately small $L$ is
essentially impossible, since the symmetry group $\mathbb{Z}_{2}^{L}$
possesses many quartic invariants. However, since we are only interested in
symmetric fixed points which are stable (viz. those which have no relevant
quartic terms), the situation becomes much more tractable. This is because we
may take advantage of powerful results about the properties of 1-loop beta
functions, which follow from interpreting the RG flow as a gradient flow on
the space of quartic couplings.

In the following we will make use of two facts. The first is that for any
subgroup $H\subset O(L)$ (we will be concerned with $H=\mathbb{Z}_{2}^{L}$),
there is at most one stable $H$-symmetric fixed point.Michel (1984); Zinn-
Justin (2002) The second is that the action of $O(L)$ on the $g_{i,j}$
couplings maps fixed points to fixed points, since for the purposes of
calculating the beta functions $O(L)$ transformations are simply redundant
relabelings of the fields.

These two facts mean that if the $g_{i,j}$ couplings describe a stable
symmetric fixed point, the $g_{i,j}$ must be invariant under the action of any
$O(L)$ transformation which preserves the $\mathbb{Z}_{2}^{L}$ symmetry
(otherwise the fixed point would not be unique). In particular, consider the
action of $\sigma\in S_{L}\subset O(L)$ on a given
$\mathbb{Z}_{2}^{L}$-invariant coupling $g_{i,j}$. The action of $\sigma$
takes $g_{i,j}\mapsto g_{\sigma(i),\sigma(j)}$, and therefore a necessary
condition for $g_{i,j}$ to give a $\mathbb{Z}_{2}^{L}$-stable fixed point is
for $g_{\sigma(i),\sigma(j)}=g_{i,j}$ for all $\sigma\in S_{L}$. This means
that any symmetric stable fixed point will have couplings of the form

$g_{i,j}^{*}=(g_{0}-h)\delta_{i,j}+h$ (103)

for some constants $g_{0},h$. If $h\neq 0$ this set of couplings gives very
non-local interactions, as it contains all-to-all quartic interlayer
couplings. On physical grounds we may then be justified in restricting our
attention to $h=0$. However, even in the general case with nonzero $h$, we
will see that no choice of $g_{0},h$ gives a stable fixed point.

Indeed, it is not hard to explicitly compute the one-loop $\beta$ functions at
fixed points of the form (103) and show that as long as $L>4$, none of them
are stable (at least within the context of the $\varepsilon$ expansion). We
find

$\displaystyle\beta_{g_{0}}$
$\displaystyle=g_{0}-\frac{9}{2}g_{0}^{2}-\frac{L-1}{2}h^{2}$ (104)
$\displaystyle\beta_{h}$ $\displaystyle=h-3g_{0}h-\frac{L+2}{2}h^{2}.$

There are three fixed points to these equations. One is the decoupled fixed
point where $h=0$, which is unstable. The other two are the $O(N)$ symmetric
fixed point $\mathcal{S}$ and the cubic fixed point $\mathcal{C}$. The
couplings at each are given by

$(g_{0}^{*},h^{*})=\begin{dcases}\mathcal{C}\,:\,&\left(\frac{2}{9}(1-1/L),\,\frac{2}{3L}\right)\\\
\mathcal{S}\,:\,&\left(\frac{2}{8+L},\,\frac{2}{8+L}\right)\end{dcases}$ (105)

It is then straightforward to check that both fixed points are unstable
(provided that $L>4$).

The above discussion has focused on the RG flows obtained in a perturbative
expansion about the free fixed point. One might imagine a possible way out by
first introducing a strong deformation that drives the system to a different
fixed point, around which the RG analysis is modified. The simplest
possibility is to add the term

$\delta S=\sum_{i}g_{pair}\int
d^{2}x\,d\tau\,\varepsilon_{2i}\varepsilon_{2i+1}.$ (106)

In the absence of further couplings between pairs of layers, this term drives
the system to a decoupled stack of 2d XY models.151515Indeed, two 2d Ising
CFTs coupled by their energy operators flow in the IR to the 2d XY fixed
point, which can be derived either from the $\varepsilon$ expansion or from
conformal perturbation theory, using the known valuesPoland et al. (2019) for
the OPE coefficients at the Ising fixed point.

This fixed point is however also unstable. Indeed, as was mentioned in the
main text, the energy operator
$\varepsilon_{i,XY}\sim\varepsilon_{2i}+\varepsilon_{2i+1}$ in the 2d XY model
is known to have scaling dimension $\Delta_{\varepsilon}\approx
1.51>3/2$,Poland et al. (2019) and as such energy-energy couplings between the
XY layers are irrelevant. However, the $\mathbb{Z}_{2}^{L}$-invariant
operators $t_{i}^{XY}=\varepsilon_{2i}-\varepsilon_{2i+1}$ (which in a given
XY layer is $\sim\phi^{2}+(\phi^{*})^{2}$) has dimension $\Delta_{t}\approx
1.24<3/2$,Chester et al. (2020) and hence couplings between $t_{i}$ operators
on neighboring layers are relevant, destabilizing the fixed point.

Therefore while we do not have a rigorous proof, the phase transition that
occurs in the $N=2$ case generically seems to be driven first-order by
fluctuations, as was suggested by the numerics of Ref. Slagle and Kim, 2017a.

##  Appendix C Effects of the ring-exchange term

In this appendix we discuss what happens to the critical point (80) in the
limit where the chains are decoupled and the ring-exchange term is the
dominant relevant perturbation (as discussed in section V.3, for $N>4$ there
is always a region where the ring-exchange term is the only relevant
perturbation to the decoupled fixed point).

When the ring-exchange term is relevant, we are prompted to expand the cosine
as $1-\frac{1}{2}a^{2}(\partial_{x}\partial_{y}\Phi)^{2}$, where now the field
$\Phi$ fluctuates in all four spacetime directions. Doing this and integrating
out the $\Theta_{w}$ fields161616The vertex operators $e^{i\Theta_{w}}$ create