# Weak symmetry breaking and topological order in a 3D compressible quantum
liquid

Joseph Sullivan  Department of Physics, Yale University, New Haven, CT
06511-8499, USA  Arpit Dua  Department of Physics and Institute for Quantum
Information and Matter, California Institute of Technology, Pasadena,
California 91125, USA  Meng Cheng  Department of Physics, Yale University, New
Haven, CT 06511-8499, USA

###### Abstract

We introduce a new type of 3D compressible quantum phase, in which the U(1)
charge conservation symmetry is weakly broken by a rigid string-like order
parameter, and no local order parameter exists. We show that this gapless
phase is completely stable and described at low energy by an infinite-
component Chern-Simons-Maxwell theory. We determine the emergent symmetry
group, which contains U(1) 0-form planar symmetries and an unusual subgroup of
the dual U(1) 1-form symmetry supported on cylindrical surfaces. Through the
associated ’t Hooft anomaly, we examine how the filling condition is fulfilled
in the low-energy theory. We also demonstrate that the phase exhibits a kind
of fractonic topological order, signified by extensively many different types
of topologically nontrivial quasiparticles formed out of vortices of the weak
superfluid. A microscopic model realizing the weak superfluid phase is
constructed using an array of strongly coupled Luttinger liquid wires, and the
connection to the field theory is established through boson-vortex duality.

###### Contents

  1. I Introduction
  2. II Coupled wire construction
    1. II.1 Model Hamiltonian and symmetries
    2. II.2 Spectrum and excitations
    3. II.3 String superfluid order
    4. II.4 Duality mapping
      1. II.4.1 OBC in z direction
    5. II.5 Vortices and quasiparticles
  3. III Infinite-component CS theory
    1. III.1 Global symmetry and anomaly
      1. III.1.1 Global symmetry
      2. III.1.2 ’t Hooft anomaly
    2. III.2 The filling anomaly
    3. III.3 Monopole operators
    4. III.4 Classification of charge excitations
  4. IV Summary and discussion
  5. A String operators for solitonic excitations
  6. B Details of the duality mapping
    1. B.1 Boson-vortex duality
    2. B.2 Integrating out matter fields
    3. B.3 Coupling to vortices
  7. C Quantization condition for 1-form symmetry transformation
  8. D Phenomenologies
    1. D.1 Order parameter
    2. D.2 Compressibility
    3. D.3 Electromagnetic response
    4. D.4 Interaction between vortices

##  I Introduction

A quantum state is compressible if the charge density can be tuned
continuously by changing the chemical potential in a translation-invariant
system. Common examples of stable compressible matter include:

  1. 1.

Symmetry-breaking phase, in particular when the U(1) charge conservation is
spontaneously broken by a local order parameter, i.e. a superfluid.

  2. 2.

For fermionic systems, the Fermi liquid is a compressible phase which does not
break any symmetry (U(1) or translation).

Both can be “deformed” to obtain other compressible phases with distinct low-
energy dynamics. For example, a Fermi liquid can be coupled to gapless bosonic
excitations (e.g. critical fluctuations or gauge fields), leading to various
kinds of non-Fermi liquids. For a superfluid, by adding additional
interactions one can change the dispersion of the Goldstone boson to alter the
low-energy physics. Examples of this include the quantum Lifshitz liquid [1]
and Bose-Luttinger liquid [2, 3].

While it is a formidable task to generally classify the low-energy dynamics of
gapless phases, one can first characterize them at the kinematic level, in
particular their emergent symmetries and the associated quantum anomalies.
Recently Ref. [4] has taken this approach to study the low-energy theory of
compressible phases. The filling condition can be formulated in a way similar
to the ’t Hooft anomaly [5, 6], and in particular it is shown that the
emergent symmetry group cannot be a compact 0-form and/or finite higher-form
symmetry group. The previous two classes of examples illustrate two mechanisms
to satisfy the filling condition: the symmetry group of the superfluid phase
is $\mathrm{U}(1)^{[0]}\times\mathrm{U}(1)^{[D-1]}$ [7, 8], where $D$ is the
spatial dimension. The presence of the $\mathrm{U}(1)^{[D-1]}$ form symmetry
is essentially equivalent to the fact that superfluid vortices must form
$(D-1)$-dimensional closed surface. The mixed anomaly between the two
subgroups is responsible for the filling condition. The same mechanism
underlies some more exotic examples of compressible phases of bosons, such as
the Bose-Luttinger liquid [3]. For the Fermi liquid, the emergent symmetry is
the (0-form) loop group[4] $\mathrm{L}^{D-1}\mathrm{U}(1)$, which describes at
the level of kinematics the existence of a Fermi surface. The group is
“larger” than any compact Lie group. Deformations of the two types of examples
may change the emergent symmetry group, but do not affect the way that the
filling anomaly is realized. In a sense, these two examples represent two
“minimal” emergent symmetry groups that allow compressibility in the
respective classes. It is clearly an important problem to understand whether
there exist other low-energy kinematical structures compatible with a generic
filling factor.

Ref. [4] mostly operates within the framework of quantum field theory in
continuum with both translation and rotation symmetries. The aforementioned
examples can be easily modified to introduce certain anisotropies, for example
in boson velocity or the shape of the Fermi surface. However these
modifications can all be smoothly turned off without changing the nature of
the states i.e. no phase transitions. In particular, they do not affect the
kinematical structure. A natural question then is whether there are
fundamentally new types of symmetry mechanisms to allow arbitrary filling,
once we go beyond the framework of continuum field theory. In fact, there is a
trivial toy example that goes beyond what we have discussed so far: consider a
stack of decoupled 2D superfluids 111Similarly, one can consider a stack of 2D
Fermi liquids. However, such a state has an open Fermi surface consisting of
two curves wrapping around the Brillouin zone.. Such a phase is obviously
compressible but the emergent symmetry group is very different from that of a
3D superfluid. More specifically, each layer has a separate $\mathrm{U}(1)$
0-form symmetry and $\mathrm{U}(1)$ 1-form symmetry at low energy. One may
object that a stack of 2D superfluids is unstable to infinitesimal boson
tunneling, thus does not represent a truly stable phase. Quantum phases with
such “subsystem” symmetries have been vigorously studied in recent years,
largely inspired by the discovery of fracton topological order [10, 11, 12,
13, 14, 15, 16, 17]. It has been recognized recently that a large class of
two- and three-dimensional quantum phases cannot be described as conventional
quantum field theories [18, 19, 20, 21, 22, 23, 24, 25, 26, 27]. They are
often characterized by exotic global symmetries [15, 28, 29], e.g. those
defined on certain submanifolds.

Motivated by this question and building on top of recent works on fracton
models [30, 31, 32], in this work we study a new example of a completely
stable compressible state in three dimensions, whose low-energy sector is
described by an infinite-component Chern-Simons (iCS) theory, consisting of
coupled planar U(1) gauge fields (other generalizations of CS-like theory to
3D in the context of fractonic order have been considered in [33] and [34]).
Variants of the field theory with fully gapped spectrum have been studied in
Ref. [30], which are shown to possess a particular kind of fracton topological
order, with all quasiparticles restricted to move in planes. Moreover, it is
proposed that the theory provides an example of a “non-foliated” fracton
topological phase. Ref. [30] also noticed that the field theory we obtain has
a gapless spectrum. One motivation for this work is to elucidate the rich
physics of the gapless iCS theory, and provide a lattice realization which can
be controllably solved and which manifests the symmetries of the underlying
field theory.

Below we briefly summarize the main features of this compressible quantum
liquid, focusing on those that distinguish it from the known types of
compressible phases previously mentioned. Most importantly, the emergent
symmetry group (relevant for the filling anomaly) consists of a 0-form
$\mathrm{U}(1)$ symmetry group of charge conservation, and a “cylindrical”
1-form symmetry, the meaning of which will be defined later. (Note that an
ordinary 3D superfluid has an emergent U(1) 2-form symmetry.) This exotic
emergent 1-form symmetry stems from a rigid string-like superfluid order,
where an order parameter for the U(1) 0-form symmetry is supported on a
straight line penetrating the entire system. However, there exists no other
local order parameter. Thus we dub this phenomenon “weak symmetry breaking”
(WSB). The terminology has been previously used in the context of (2+1)d
symmetry-enriched topological order [35, 36, 37], when there exists no local
order parameter but certain non-local observables break the (0-form) symmetry.
In addition, the low-energy theory also exhibits finite 1-form symmetries,
corresponding to deconfined, topologically nontrivial quasiparticle
excitations. They can be formed from dipoles of vortices of the “weak”
superfluid. Thus in a sense the phase of matter intertwines U(1) symmetry
breaking and topological order in an interesting way.

Phenomenologically, we find that at low energy there is only a single gapless
point. We argue that the gapless state is in fact robust to any local
perturbation and therefore represents a stable gapless phase. In particular,
the gapless excitations are charge neutral and all local charged excitations
are gapped. We also study the transport properties, and find that the phase
exhibits a superconducting response in the plane perpendicular to the
direction of the non-local order parameter, and insulating in the other.

The paper is organized as follows: in Sec. II we first provide a microscopic
realization of the compressible liquid with WSB, using a coupled-wire
construction, and study its various properties. We derive an infinite-
component CS theory description of the low-energy physics, by applying boson-
vortex duality transformation in the coupled wire setting. Gapped cases of
these field theories have been recently studied in Ref. [30] as examples of
Abelian fractonic phases without any foliation structure. In Sec. III we then
explore the emergent symmetries, anomalies (especially the filling anomaly)
and stability of the phase, using the dual gauge theory formulation.

##  II Coupled wire construction

We start from a microscopic model realizing the phase. The model is built from
an array of interacting one-dimensional quantum wires. Such models have proven
very fruitful in providing explicit realizations of topological phases in both
two and three dimensions [38, 39, 31, 32, 40, 32, 31].

###  II.1 Model Hamiltonian and symmetries

Consider quantum wires arranged in a square lattice. Each wire is described by
a $c=1$ bosonic Luttinger liquid, with a K matrix $K_{\text{w}}=\sigma^{x}$.
We postulate that the Luttinger liquid is realized as the low-energy effective
theory of a one-dimensional chain of bosons (or spins). The Hamiltonian is

$H=\frac{v}{2\pi}\sum_{\mathbf{r}}\int\mathrm{d}x\,[(\partial_{x}\varphi_{\mathbf{r}})^{2}+(\partial_{x}\theta_{\mathbf{r}})^{2}].$
(1)

Here ${\mathbf{r}}=(y,z)$ labels the position of wires in the $yz$ plane. The
bosonic fields $\varphi$ and $\theta$ satisfy the canonical commutation
relation

$[\varphi_{\mathbf{r}}(x_{1}),\partial_{x_{2}}\theta_{\mathbf{r}^{\prime}}(x_{2})]=2\pi
i\delta(x_{1}-x_{2})\delta_{\mathbf{r}\mathbf{r}^{\prime}}.$ (2)

We have also assumed that all wires have the same velocity to ensure
translation invariance.

We add the following type of interactions to gap out the wires:

$-g\sum_{\mathbf{r}}\int\mathrm{d}x\,\cos\Theta_{\mathbf{r}}(x),g>0.$ (3)

Here $\Theta_{\mathbf{r}}\equiv\Theta_{y,z}$ is defined as

$\begin{split}\Theta_{y,z}&=-\varphi_{yz}+m\theta_{yz}+\varphi_{y+1,z}+m\theta_{y+1,z}\\\
&+(n_{1}\theta_{y,z-1}+n_{2}\theta_{y,z+1}+n_{2}\theta_{y+1,z-1}+n_{1}\theta_{y+1,z+1}).\end{split}$
(4)

This is a generalization of the coupled wire construction for bilayer quantum
Hall states by Teo and Kane [39]. The interaction $\cos\Theta$ is illustrated
in Fig. 1(b). One can easily show that these fields satisfy the null vector
condition[41]:

$[\Theta_{\mathbf{r}}(x),\Theta_{\mathbf{r}^{\prime}}(x^{\prime})]=0,$ (5)

so they can be minimized simultaneously. We are interested in the strong-
coupling limit $g\rightarrow\infty$, at which the cosine terms pin all
$\Theta$ fields to the minima. This can be achieved either by having a large
bare value of $g$, or turning on inter-wire density-density interactions to
make $g$ a relevant coupling.

Figure 1: The coupled wire construction in this work starts from a 2D array of
single-component Luttinger liquids, as illustrated in (a), and then gapping
interactions between nearby wires are turned on to create a 3D quantum state.
The form of the interaction is illustrated in (b).

First we study the global symmetry of the model. Before turning the coupling
Eq. (4), each wire has $\mathrm{U}(1)_{\varphi}\times\mathrm{U}(1)_{\theta}$
symmetry, the charge densities of which are $\frac{1}{2\pi}\partial_{x}\theta$
and $\frac{1}{2\pi}\partial_{x}\varphi$, respectively, and there is a mixed ’t
Hooft anomaly between the two $\mathrm{U}(1)$ subgroups. One of them,
$\mathrm{U}(1)_{\varphi}$, can be understood as the conservation of boson
number, which is assumed to be an exact symmetry of the microscopic
Hamiltonian. Here we choose the convention that the boson density is given by
$\frac{1}{2\pi}\partial_{x}\theta$, and the symmetry transformation acts on
the fields as

$\varphi\rightarrow\varphi+\alpha,$ (6)

where $\alpha\in[0,2\pi)$ is the $\mathrm{U}(1)$ rotation angle. The other
subgroup $\mathrm{U}(1)_{\theta}$ then must be emergent at low energy due to
the mixed anomaly. Importantly, the lattice translation symmetry along the
wire is realized at low energy as a particular element of
$\mathrm{U}(1)_{\theta}$:

$T_{x}:\theta\rightarrow\theta+2\pi\nu,$ (7)

where $\nu$ is the filling factor of the wire, which is equal to the filling
factor of the entire 3D system. When $\nu$ is not an integer, the celebrated
Lieb-Schultz-Mattis theorem [42] shows that the 1D wire cannot be fully gapped
without breaking the translation symmetry.

The coupling Eq. (4) apparently reduces the symmetry group. For each $xy$
plane, there is a $\mathrm{U}(1)_{\varphi}$ subsystem symmetry:

$\varphi_{yz}\rightarrow\varphi_{yz}+\alpha_{z},$ (8)

However, the $\mathrm{U}(1)_{\theta}$ symmetry is broken in general.

We are mainly interested in the special case when $m=-(n_{1}+n_{2})$. For this
choice of parameters, we find a $\mathrm{U}(1)_{\theta}$ subsystem symmetry in
each $xz$ plane:

$\theta_{yz}\rightarrow\theta_{yz}+\beta_{y}.$ (9)

This obviously includes a global transformation
$\theta_{\mathbf{r}}\rightarrow\theta_{\mathbf{r}}+\beta$.

If we set $n_{1}=n_{2}$, then there is a further U(1) dipole symmetry on each
$xz$ plane:

$\theta_{yz}\rightarrow\theta_{yz}+z\gamma_{y}.$ (10)

We now consider spatial translations. The full Hamiltonian is obviously
invariant under the discrete translations along $y$ and $z$ directions. For
the lattice translation $T_{x}$, the interaction is invariant if and only if
$(m+n_{1}+n_{2})\nu$ is an integer. If $m=-(n_{1}+n_{2})$, then $\nu$ can be
any real number. Otherwise the model is only translation-invariant for certain
special rational fillings. We thus conclude that the model with
$m=-(n_{1}+n_{2})$ can be defined _at any filling_. As will be shown below,
there is no spontaneous symmetry breaking for either the global
$\mathrm{U}(1)_{\varphi}$ or $\mathrm{U}(1)_{\theta}$, in the sense that no
local order parameter exists. Therefore the coupled wire model constitutes a
new example of compressible quantum phase.

To summarize, we find that at $m=-(n_{1}+n_{2})$, the system preserves the
following global continuous symmetry:

$\begin{gathered}\mathrm{U}(1)_{\varphi}:\varphi_{\mathbf{r}}\rightarrow\varphi_{\mathbf{r}}+\alpha,\\\
\mathrm{U}(1)_{\theta}:\theta_{\mathbf{r}}\rightarrow\theta_{\mathbf{r}}+\beta.\end{gathered}$
(11)

These two U(1) symmetries have a mixed ’t Hooft anomaly, so cannot be realized
as on-site symmetries microscopically at the same time. In this work we choose
the convention that $\mathrm{U}(1)_{\varphi}$ is the boson number
conservation, and $\mathrm{U}(1)_{\theta}$ is an emergent symmetry. The
lattice translation along wires is embedded into $\mathrm{U}(1)_{\theta}$.

So far we have only considered the so-called 0-form symmetries. At low energy,
the system can develop emergent “higher-form” symmetries [43], whose charges
are extended objects. An example relevant to our discussion is the U(1) 1-form
symmetry in a 2D superfluid, which is physically equivalent to all superfluid
vortices being non-dynamical (i.e. infinitely heavy). In our model, if
$m=n_{1}=n_{2}=0$, the system is a stack of layers of 2D superfluids, and at
low energy each layer has its own U(1) 1-form symmetry. Schematically, the
1-form charge in the layer $z$ for a closed path $C$ is given by

$Q_{z}(C)=\oint_{C}\nabla\varphi_{z},$ (12)

which counts the total vorticity enclosed by $C$. When the coupling to
$\theta$ is turned on, since $e^{i\theta}$ creates a vortex-anti vortex pair
and hops vortices in the layer, the U(1) 1-form symmetry for the layer is
explicitly broken. However, when $m=-n_{1}-n_{2}$, the term
$\cos(\Theta_{y,z})$ preserves the 1-form charge

$Q_{z-1}(C)+Q_{z}(C)+Q_{z+1}(C).$ (13)

Note that $C$ must be the same for the three layers. Now considering all
layers together, the total 1-form charge

$\sum_{z}Q_{z}(C)$ (14)

is conserved by the Hamiltonian at low energy when the cosine term dominates,
assuming periodic boundary condition. We refer to the conservation of the
total vorticity Eq. (14) as a cylindrical 1-form symmetry, since the closed
loop $C$ must be exactly aligned throughout all the layers, so the symmetry
operator is in fact supported on a cylindrical surface. This should be
compared with the 2-form symmetry in a 3D superfluid, where the symmetry
operator is supported on arbitrary loops. We illustrate the symmetry operator
in Fig 2b.

Figure 2: (a) The 2-form symmetry operators in a 3D superfluid are supported
on closed curves. The corresponding charged objects form closed surfaces (i.e.
vortex sheets), hence a 2-form symmetry. (b) The higher-form symmetry operator
of the WSB compressible liquid, which is defined in Eq 14. The surface is
deformable in $xy$-plane provided it is deformed in the same way in each
layer. We refer to this as a cylindrical 1-form symmetry.

###  II.2 Spectrum and excitations

We now move on to analyze the spectrum of the Hamiltonian. The model is
generally not solvable, so we limit ourselves to the strong-coupling limit
where the cosine potential terms dominate over the kinetic terms. As shown in
Ref. [31], low-energy excitations can be divided into two types: first of all,
there are excitations that describe smooth fluctuations of the bosonic fields,
similar to spin waves. They can be analyzed in the “mean-field” approximation,
where we expand $\cos\Theta_{\mathbf{r}}\approx
1-\frac{1}{2}\Theta_{\mathbf{r}}^{2}$, and then solve the quadratic theory. We
find that the spectrum of such Gaussian fluctuations is given by

$E_{\mathbf{k}}=\sqrt{v^{2}k_{x}^{2}+vg(f_{\varphi}+f_{\theta})}$ (15)

where

$\begin{split}f_{\varphi}&=2-2\cos k_{y},\\\
f_{\theta}&=4\left[m\cos\frac{k_{y}}{2}+n_{1}\cos(\frac{k_{y}}{2}-k_{z})+n_{2}\cos(\frac{k_{y}}{2}+k_{z})\right]^{2}.\end{split}$
(16)

When $|m|\leq|n_{1}+n_{2}|$, one can make both $f_{\varphi}=f_{\theta}=0$ by
setting $k_{y}=0,\cos k_{z}=-\frac{m}{n_{1}+n_{2}}$, which is clearly the
minimum of $f_{\varphi}+f_{\theta}$ since both functions are non-negative.
Therefore the spectrum is gapless. For $m=-(n_{1}+n_{2})$, the minimum is at
$\mathbf{k}=(0,0,0)$. Near this point, $E_{\mathbf{k}}$ takes the following
approximate form:

$\begin{split}E_{\mathbf{k}}&=\Big{[}v^{2}k_{x}^{2}+vg\Big{(}k_{y}^{2}+2(n_{2}^{2}-n_{1}^{2})k_{y}k_{z}^{3}\\\
&\quad\quad+(n_{1}-n_{2})^{2}k_{y}^{2}k_{z}^{2}+(n_{1}+n_{2})^{2}k_{z}^{4}\Big{)}\Big{]}^{1/2}.\end{split}$
(17)

For $n_{1}=n_{2}=n$, we have

$f_{\varphi}+f_{\theta}=4\cos^{2}\frac{k_{y}}{2}[(m+2n\cos k_{z})^{2}-1]+4.$
(18)

If $|m|>2|n|$, the minimum of this function is 4 at $k_{y}=\pm\pi$ and the
spectrum is fully gapped. When $|m|\leq 2|n|$, the minimum is $0$ with
$k_{y}=0$ and $\cos k_{z}=-\frac{m}{2n}$. In this case the spectrum is
gapless, with one or two gap-closing points. For $m=-2n$, near the gap-closing
point $\mathbf{k}=(0,0,0)$, the spectrum is approximately

$E_{\mathbf{k}}=\sqrt{v^{2}k_{x}^{2}+vgk_{y}^{2}+4vgn^{2}k_{z}^{4}}.$ (19)

Besides smooth fluctuations, there also exist localized excitations,
corresponding to “discontinuities” in the field configurations. At the mean-
field level, they can be built out of soliton excitations of the gapping terms
that tunnel between different minima. More concretely, for
$\cos\Theta_{\mathbf{r}}(x)$, a $k$-soliton where $k\in\mathbb{Z}$ at $x_{0}$
is a configuration where $\Theta_{\mathbf{r}}$ winds by $2\pi k$ over a short
distance $\xi$ around $x_{0}$. In this limit, these excitations can be thought
of as massive quasi-particles. While we do not know the exact profile of such
an excitation, we assume that they stay gapped beyond the mean-field analysis,
at least for a range of $g$. We give a more systematic description of the
solitonic excitations in Sec. II.5.

To create such excitations, one can apply vertex operators, such as
$e^{i\varphi}$ or $e^{i\theta}$. Note that these two, together with
derivatives of $\varphi$ and $\theta$, generate all local operators in the
wire model. It is easy to show that derivatives of $\varphi$ or $\theta$
cannot create localized solitons. Rather they only result in smooth
configurations of the fields. Therefore we can focus just on the vertex
operators. A crucial fact which we rely on here is that the model satisfies a
version of the “topological order” condition as discussed in Ref. [31].
Namely, if a vertex operator has a finite support and commutes with all
$\Theta_{\mathbf{r}}$, then it must be a linear combination of
$\Theta_{\mathbf{r}}$’s. In other words, any vertex operator with local
support has to create gapped excitations, provided it is not some combination
of $\Theta_{\mathbf{r}}$’s. Therefore such operators are infinitely irrelevant
even when the theory is gapless. This strongly suggests that the gapless
theories, when $|m|\leq|n_{1}+n_{2}|$, are stable with respect to weak local
perturbations, as long as the topological order condition is obeyed. This is
analogous to 3D U(1) gauge theory with gapped charges.

It is important to note that the topological order condition is shown for an
infinitely large system. When the system is finite (or finite in one
direction), the situation turns out to be quite different. As we discuss next,
there can be vertex operators commuting with all gapping terms. Such operators
are crucial in understanding the symmetry breaking (or the absence thereof) in
the low-energy theory.

###  II.3 String superfluid order

In this section we focus on $m=-(n_{1}+n_{2})$ and assume that both global
$\mathrm{U}(1)_{\theta}$ and $\mathrm{U}(1)_{\varphi}$ symmetries are
preserved. The topological order condition does not preclude the existence of
non-local vertex operators that commute with all plaquette terms. We can
consider a closed system with periodic boundary condition along $z$, in which
case we can identify two such non-local operators:

$\Phi_{y}=\sum_{z}\varphi_{yz},~{}\Pi_{z}=\sum_{y}\theta_{yz}.$ (20)

Note that when $n_{1}=n_{2}$ there is an additional operator of interest given
by $D_{y}=\sum_{z}z\varphi_{yz}$, but it does not exist with PBC. The operator
$e^{i\Phi_{y}}$ has a charge $N_{z}$ under $\mathrm{U}(1)_{\varphi}$, where
$N_{z}$ is the size of system in the $z$ direction. Likewise, $e^{i\Pi_{z}}$
has a charge $N_{y}$ under $\mathrm{U}(1)_{\theta}$.

We show in Appendix D.1 that with a finite $N_{z}$, when viewed as a quasi-
two-dimensional system, the $e^{i\Phi_{y}}$ operator orders and spontaneously
breaks the $\mathrm{U}(1)_{\varphi}$ symmetry. More precisely, the two-point
function decays as

$\langle e^{i\Phi_{y}(x)}e^{-i\Phi_{y^{\prime}}(x^{\prime})}\rangle\rightarrow
e^{-\frac{L_{z}}{\rho|\mathbf{r}-\mathbf{r}^{\prime}|}},$ (21)

where $\rho$ is a constant, and $\mathbf{r}=(x,y)$ is the coordinate in the xy
plane. Therefore, fixing $N_{z}$, the correlation function approaches a
constant when $|\mathbf{r}-\mathbf{r}^{\prime}|$ is much greater than
$N_{z}/\rho$. This makes sense since for such large separation the thickness
along the $z$ direction becomes negligible. However, going to the 3D limit
when $N_{z}$ is comparable to the size of the other dimensions, the superfluid
order is suppressed and eventually vanishes in the thermodynamic limit
$N_{z}\rightarrow\infty$.

Equivalently, one can consider the superfluid stiffness in the quasi-two-
dimensional system. The “Josephson” coupling that stabilizes the superfluid
order originates from the $\varphi_{y,z}-\varphi_{y+1,z}$ term in
$\Theta_{\mathbf{r}}$. However, such a term only accesses a small local
portion of the order parameter, which is highly non-local and runs through the
entire $z$ direction. As a result, one can show that the superfluid stiffness
vanishes as $1/N_{z}$.

To summarize, we find that the line operator $e^{i\Phi_{y}}$ develops long-
range order breaking the U(1)φ symmetry, but there exist no local U(1)φ order
parameter. This is what we refer to as the phenomenon of “weak symmetry
breaking”. We note that similarly $e^{i\Pi_{z}}$ exhibits long-range order in
the $x-z$ plane, although the asymptotic of the correlation function is now
more anisotropic. So the U(1)θ symmetry is also weakly broken by the line
order parameter.

It is also instructive to consider open boundary conditions along $z$. Suppose
there are two $yz$ surfaces at $z=0$ and $z=N_{z}-1$. We first choose the
gapping terms to be

$-g\sum_{y}\sum_{z=1}^{N_{y}-2}\cos\Theta_{yz}.$ (22)

At each $y$, there are $N_{z}$ wires but only $N_{z}-2$ gapping terms, which
leaves certain fields near the boundary unpinned. We focus on the $z=0$
surface. Unpinned fields on the surface are generated by $\theta_{y0}$ and

$\tilde{\varphi}_{y+\frac{1}{2}}=\varphi_{y+1,0}-\varphi_{y,0}+n_{1}\theta_{y,1}+n_{2}\theta_{y+1,1}.$
(23)

Note that both $\tilde{\varphi}_{y+\frac{1}{2}}$ and $\theta_{y0}$ commute
with themselves. We see that all unpinned fields are neutral under
$\mathrm{U}(1)_{\varphi}$ symmetry.

Now we further add interactions on the surface, e.g.
$\cos(\tilde{\varphi}_{y-\frac{1}{2}}-\tilde{\varphi}_{y+\frac{1}{2}})$ or
$\cos(\theta_{y0}-\theta_{y+1,0})$, which spontaneously break the
$\mathrm{U}(1)_{\theta}$ symmetry. In this case, the $e^{i\Phi_{y}}$ operator
no longer commutes with the surface order and thus is not a low-energy
observable.

We can also create a completely symmetric and gapped surface, by adding the
following perturbations:

$-\cos\left[\frac{m}{2}(\theta_{y,0}+\theta_{y+1,0})+\tilde{\varphi}_{y+\frac{1}{2}}\right].$
(24)

One can check that all the cosine arguments commute with each other and there
are enough of them to gap out all the degrees of freedom on the surface. In
this case, one can slightly modify the definition of $\Phi_{y}$ at the
surfaces so that it still commutes with all the gapping interactions. This is
done by adding a vertex operator which generates an appropriate phase slip at
the end of the string operator:
$\Phi_{y}\to\sum_{z}\varphi_{yz}-(\frac{n_{1}-n_{2}}{2})\theta_{y0}$. Note
that when $n_{1}=n_{2}=-m/2$, $\Phi_{y}$ has the same form as the PBC case.

While we do not have a general proof, we believe one of following scenarios
must occur: either the $\mathrm{U}(1)_{\theta}$ symmetry is spontaneously or
explicitly broken at the surface, in which case no (local or non-local)
$\mathrm{U}(1)_{\varphi}$ order parameter exists in the system, or the
$\mathrm{U}(1)_{\theta}$ symmetry is preserved but one can then find a non-
local order parameter by modifying $\Phi_{y}$ near the surface.

Since the model has a non-local superfluid order, one may wonder whether there
is any superconducting response once coupled to an external electromagnetic
field. Given the model is highly anisotropic, it is not surprising that the
answer depends on direction. As shown in Eq. (8), the U(1) charge in each
layer is conserved by the Hamiltonian as well as the ground state. As a result
there cannot be any charge transport in the $z$ direction, i.e. it is
insulating. The other directions are quite different. Using the dual gauge
theory we compute the effective action for the external field in Appendix D.3,
and indeed show that the system is superconducting in the $xy$ plane.

###  II.4 Duality mapping

In this section we describe how the microscopic wire model can be mapped to an
infinite-component CS theory coupled to gapped matter fields. The full details
of this procedure are presented in Appendix B. Schematically, the gauge theory
emerges by employing a boson-vortex duality within each horizontal layer of
wires. We adopt the method explained in Refs. [44,45] to the 3D coupled wire
model. Note that an alternative method to derive Chern-Simons-type gauge
theory from coupled wire models have been proposed in Ref. [46].

Recall that the wire model is described by a rectangular array of Luttinger
liquids with conjugate variables $(\varphi,\theta)$ and Lagrangian

$\begin{split}\mathcal{L}[\varphi,\theta]&=\sum_{\mathbf{r}}\frac{i}{\pi}\partial_{x}\theta_{\mathbf{r}}\partial_{\tau}\varphi_{\mathbf{r}}+\frac{\tilde{v}}{2\pi}\left(\partial_{x}\varphi_{{\mathbf{r}}}\right)^{2}+\frac{u}{2\pi}\left(\partial_{x}\theta_{{\mathbf{r}}}\right)^{2}\\\
&~{}~{}~{}~{}~{}~{}+\frac{v}{8\pi}\left(\partial_{x}\Delta_{y}\varphi_{{\mathbf{r}}}\right)^{2}-g\cos(2\bar{\theta}_{\mathbf{r}})\end{split}$
(25)

where

$\begin{split}2\bar{\theta}_{{\mathbf{r}}}&=\Delta_{y}\varphi_{\mathbf{r}}+\Lambda\theta_{\mathbf{r}}\\\
&=\left(\varphi_{{\mathbf{r}}+\hat{y}}-\varphi_{{\mathbf{r}}}\right)\\\
&~{}~{}~{}~{}+(m\theta_{{\mathbf{r}}}+m\theta_{{\mathbf{r}}+\hat{y}}+n_{1}\theta_{{\mathbf{r}}-\hat{z}}+n_{2}\theta_{{\mathbf{r}}+\hat{z}}\\\
&~{}~{}~{}~{}+n_{2}\theta_{{\mathbf{r}}+\hat{y}-\hat{z}}+n_{1}\theta_{{\mathbf{r}}+\hat{y}+\hat{z}}).\end{split}$
(26)

Here we define $\Delta_{y}X_{\mathbf{r}}=X_{{\mathbf{r}}+\hat{y}}-X_{r}$ and
${\mathbf{r}}=(y,z)$ is a wire index. Note that the kinetic part of Eq. (25)
is more general than the free boson Hamiltonian given in Eq. (1). In
particular,
$\frac{v}{8\pi}\left(\partial_{x}\Delta_{y}\varphi_{{\mathbf{r}}}\right)^{2}$
has been added to the standard Luttinger Liquid kinetic term. This is done
because it is convenient for the duality mapping to the vortex theory but its
presence does not affect the qualitative physics [44].

Throughout the discussion of the microscopic model we use two equivalent forms
of labelling for the wires, $O_{{\mathbf{r}}+\hat{a}}\equiv
O_{y+a_{y},z+a_{z}}$. Anticipating the layered structure of the gauge theory,
we treat $z$ as the “layer” index while $y$ is coarse grained to a continuous
spatial coordinate. With this motivation in mind, we define the following pair
of conjugate variables:

$\begin{array}[]{l}\tilde{\varphi}_{yz}=-\sum_{y^{\prime}}\operatorname{sgn}\left(y^{\prime}-y-\frac{1}{2}\right)\theta_{y^{\prime}z}\\\
\\\
\tilde{\theta}_{yz}=\frac{1}{2}\left(\varphi_{y+1,z}-\varphi_{y,z}\right).\end{array}$
(27)

Here $\tilde{\varphi}_{yz}$ creates a $2\pi$ vortex in the $\varphi$ field in
layer $z$, in between wires $y$ and $y+1$. The operator
$\partial_{x}\tilde{\theta}_{yz}$ is the “charge” operator for this vortex. We
can re-express the Lagrangian in Eq. (25) in terms of these new fields but the
result will be highly non-local in the $y$ direction. We can restore locality
in each layer $z$ via a Hubbard-Stratonovich transformation at the expense of
adding new degrees of freedom $a_{0}^{(z)}(x,y)$ and $a_{1}^{(z)}(x,y)$. At
this stage, the Lagrangian takes the form

$\begin{split}\mathcal{L}[\tilde{\varphi},\tilde{\theta},a_{\mu}]&=\sum_{r}\frac{i}{\pi}\partial_{x}\tilde{\theta}_{\mathbf{r}}(\partial_{\tau}\tilde{\varphi}_{\mathbf{r}}-a_{0,{\mathbf{r}}}^{(z)})+\frac{u}{2\pi}\left(\partial_{x}\tilde{\varphi}_{{\mathbf{r}}}-a_{1,{\mathbf{r}}}^{(z)}\right)^{2}\\\
&~{}~{}~{}~{}~{}~{}+\frac{v}{2\pi}\left(\partial_{x}\tilde{\theta}_{\mathbf{r}}\right)^{2}-g\cos(2\tilde{\theta}_{\mathbf{r}}+\frac{1}{2}\Lambda\cdot\Delta_{y}\tilde{\varphi}_{\mathbf{r}})\\\
&~{}~{}~{}~{}~{}~{}+\mathcal{L}_{\text{Maxwell}}^{(z)}.\end{split}$ (28)

Here, the definition of $\Lambda$ can be easily inferred from Eq. (26), so
that the argument of the cosine term matches $2\bar{\theta}_{\mathbf{r}}$.
$\mathcal{L}_{\text{Maxwell}}^{(z)}=\frac{u}{2\pi}(\Delta_{y}a_{1}^{(z)})^{2}+\frac{1}{2\pi\tilde{v}}(\Delta_{y}a_{0}^{(z)})^{2}$.
The superscript of this term is no accident; we interpret Eq. (28) as a stack
of 2D gauge theories in the $a_{2}^{(z)}=0$ gauge coupled to matter
$(\tilde{\varphi},\tilde{\theta})$. In the context of Eq. 28 the gauge fields
would appear to be over-labeled since ${\mathbf{r}}=(y,z)$. We keep the
redundant $(z)$ superscript to emphasize this picture of stacks of 2D gauge
theories. Such “layered” gauge theories have recently been studied in the
context of fracton physics [30] and are reminiscent of the foliated gauge
theories developed to describe fracton order [47]. We emphasize that the
interpretation of $a$’s as a stack of 2D gauge fields in the $xy$ plane is
motivated, and ultimately justified by providing a correct description of the
low-energy physics of the model.

The next step is to integrate out the fields $\tilde{\varphi}$ and
$\tilde{\theta}$. The details of the derivation can be found in Appendix B. In
the end we obtain the effective layered Maxwell-CS theory:

$\mathcal{L}[a]\equiv\frac{iK_{zz^{\prime}}}{4\pi}a^{(z)}\wedge
da^{(z^{\prime})}+(\text{Maxwell terms}),$ (29)

where

$K_{zz^{\prime}}=2m\delta_{zz^{\prime}}+(n_{1}+n_{2})\delta_{z,z^{\prime}\pm
1}.$ (30)

The inter-layer Maxwell term contains couplings between components of the
electromagnetic field from adjacent layers $(z)$; as an example the $B^{2}$
term corresponds to
$\sum_{z,z^{\prime}}B^{(z)}\left[\frac{u}{2\pi}\delta_{zz^{\prime}}+\frac{v}{8\pi}(K^{\mathsf{T}}K)_{zz^{\prime}}\right]B^{(z^{\prime})}$.

So far we have assumed that $\varphi$ and $\theta$ vary smoothly, excluding
the “solitonic” excitations of the cosine pinning term. These more singular
configurations are the analog of vortices in the usual boson-vortex duality,
which should be minimally coupled to the dynamical gauge fields. In Appendix
B.3 we show that this is indeed the case; namely the solitons become gauge
charges in the dual gauge theory.

####  II.4.1 OBC in z direction

Now we consider OBC in the $z$ direction and the effect of the various surface
gapping terms discussed in Sec II.3 on the resulting dual gauge theory.

First we discuss the symmetrically gapped surface generated by the pinning
term in Eq. (24). The analysis is similar to the PBC case, except now we must
pay attention to the $z=0$ layer of wires. We decompose $\sum_{\mathbf{r}}$
into $\sum_{y}\sum_{z>0}+\sum_{y,z=0}$ in Eq (25). The terms in the first sum
are unchanged while in the second sum, $\sum_{y,z=0}$, the surface pinning
term is $\cos 2\bar{\theta}_{y0}$ where

$\begin{split}2\bar{\theta}_{y0}&=2\tilde{\theta}_{y0}+\frac{m}{4}\Delta_{y}(\bar{\varphi}_{y,0}+\bar{\varphi}_{y+1,0})\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{n_{1}}{2}\Delta_{y}\bar{\varphi}_{y,1}+\frac{n_{2}}{2}\Delta_{y}\bar{\varphi}_{y+1,1}~{}~{}~{}.\end{split}$
(31)

From here one can integrate out the matter fields and truncate higher order
derivative terms where appropriate. The upshot of this is a modification of
the K-matrix:

$\mathbf{K}_{\text{sym,OBC}}=m\begin{pmatrix}-1&1&&&0\\\ 1&-2&1&&\\\
&\ddots&\ddots&\ddots&\\\ &&1&-2&1\\\ 0&&&1&-1\end{pmatrix}.$ (32)

The $B^{2}$ term is also modified near the surface. This can be seen by noting
that the term in the Lagrangian corresponds to something proportional to
$(1+\alpha K^{\mathsf{T}}K)_{zz^{\prime}}B^{(z)}B^{(z^{\prime})}$ for some
constant $\alpha$.

Alternatively, we can explicitly break the U(1)θ symmetry by adding the very
natural pinning term corresponding to the middle and top layer of the bulk
plaquette term. This partial plaquette is given by

$\begin{split}2\bar{\theta}_{y0}&=2\tilde{\theta}_{y0}+\frac{m}{2}\Delta_{y}(\bar{\varphi}_{y,0}+\bar{\varphi}_{y+1,0})\\\
&~{}~{}~{}+\frac{n_{1}}{2}\Delta_{y}\bar{\varphi}_{y,1}+\frac{n_{2}}{2}\Delta_{y}\bar{\varphi}_{y+1,1}.\end{split}$
(33)

The duality mapping of this surface theory proceeds in the same way as the
symmetric gapped example just discussed. In this case

$\mathbf{K}_{\text{broken,OBC}}=m\begin{pmatrix}-2&1&&&0\\\ 1&-2&1&&\\\
&\ddots&\ddots&\ddots&\\\ &&1&-2&1\\\ 0&&&1&-2\end{pmatrix}.$ (34)

###  II.5 Vortices and quasiparticles

Throughout our derivation of the dual gauge theory, we have assumed that the
bosonic fields $\varphi$ and $\theta$ vary smoothly in spacetime. As we have
already mentioned, the complete spectrum also contains solitonic
configurations. We show in Appendix B.3 that the solitons are in fact
dynamical charges of the dual gauge fields. Below we describe the universal
properties of the soliton excitations. Many details are delegated to Appendix
A.

Physically, a “fundamental” soliton, i.e. a $2\pi$ jump in a single
$\Theta_{\mathbf{r}}$ term, is in fact a vortex of the weak superfluid. This
can be demonstrated by explicitly constructing a string operator to create a
pair of such solitons separated in the $y$ direction. The string operator
takes the following form:

$\prod_{y}e^{i\theta_{yz}(x)},$ (35)

which causes a branch cut of $2\pi$ phase jump in the string order parameter,
thus the corresponding excitation is a vortex. It is useful to divide the
excitations into vortices and non-vortices. The simplest non-vortex excitation
is just a dipole of two vortices, with opposite vorticities $k$ and $-k$,
separated along the $z$ direction. Such non-vortex excitations can be
considered as being “deconfined” quasiparticles in this phase. We call $k$ the
strength of the dipole.

Now we discuss the mobility of excitations. The string operator Eq. (35) shows
that all excitations can move freely along the $y$ direction. In Appendix A we
show that such a fundamental soliton is immobile in the $z$ direction, when
$|m|\geq 2|n|$. Notice that these results hold true in both gapped and gapless
models. For the gapless theory with $m=-2n$, one can show that a dipole whose
strength is a multiple of $m$ can move in the $z$ direction.

The question of mobility along the $x$ direction is more subtle. We present a
construction of a $x$-string operator in Appendix. The string operator is
ultra-local in the $y$ direction, but can have some extension in the $z$
direction. For a gapped model, we find that the construction does yield an
exponentially localized string operator, as expected. However, the gapless
case is quite different. To be concrete let us consider a finite $N_{z}$
system. For a vortex, the construction fails if PBC or the symmetric OBC is
imposed in the $z$ direction. Only when the symmetry-breaking OBC is present
one can find a string operator. For a dipole, a string operator always exists
regardless of the boundary condition in the $z$ direction. In both cases,
however, the constructed string operator decays only algebraically in the $z$
direction away from the localization of the excitation. This observation
suggests that the mobility of these excitations along the $x$ direction is
also reduced.

In a gapped system, one can define the superselection sectors of quasiparticle
excitations as the equivalence classes under local operations. Namely, two
excitations are equivalent if they can be transformed to each other by local
operator. We can attempt to generalize the notion to our model. We focus on
the non-vortex excitations, because there are infinitely many types labeled by
the vorticity. Since every non-vortex excitation is mobile along $y$, we can
just consider them at a fixed $y$. Equivalence between excitations located at
different $x$ is less clear since the string operator is only algebraically
localized. If however we consider them to be equivalent, we find that there
are $m^{N_{z}-1}\cdot N_{z}$ classes of non-vortex excitations when PBC is
imposed, which agrees with the topological ground state degeneracy of the
coupled wire model. More details of the counting can be found in Appendix A.

Since the vortices are coupled to the gapless gauge fields, gauge fluctuations
induce long-range interactions between the vortices. We compute the
interactions when the simplest, most isotropic Maxwell term is added to the
dual CS theory (see Eq. (57) below). For two unit vortices separated by
$\mathbf{r}=(x,y,z)$, let $\rho=\sqrt{x^{2}+y^{2}}$. The interaction potential
takes the following asymptotic form:

$V(\mathbf{r})\sim\begin{cases}\frac{1}{|z|}&\tilde{g}z^{2}\gg\rho\\\
\frac{1}{\sqrt{\rho}}&\tilde{g}z^{2}\ll\rho\end{cases}.$ (36)

From this one can further obtain the interactions between dipoles, which decay
with a higher power. Details of the calculations can be found in Appendix D.4.

##  III Infinite-component CS theory

Our duality mapping suggests that the low-energy physics can be described by
an infinite-component Chern-Simons gauge theory, as given in Eq. (29). The
gapped case $|m|>|n_{1}+n_{2}|$ was studied in greater detail in Ref. [30].
Briefly, such a theory exhibits a kind of fracton topological order, with all
quasiparticles being planons. Interestingly, in many cases the Wilson loop
operators for these quasiparticles must have exponentially decaying tails in
the $z$ direction, and as a result, the braiding statistics between
quasiparticles are not strictly local. When this phenomenon happens, the
fractonic order appears to be beyond the topological defect network framework
[48, 49, 50].

In the following we focus more on the gapless case, and our main task is to
analyze the symmetry, anomaly and the issue of stability. We show that the iCS
theory indeed captures the universal aspects of the coupled wire model. First
of all, we note that the photon spectrum of the iCS theory with only “intra-
layer” Maxwell terms was calculated in Ref. [30] and the result agrees 222At
first glance the spectrum of the layered gauge theory arising from Eq. (87)
would seem to depend only on $n_{1}+n_{2}$ and $m$, which does not agree with
the spectrum derived for the wire model. This discrepancy stems from the fact
the $y$ direction has been coarse grained in the layered picture. Restoring
higher order derivative terms in $y$ in the gauge theory gives a spectrum
matching Eq. (17)  with the Gaussian spectrum Eq. (17) of the coupled wire
model. Below we discuss how other universal aspects of the coupled wire model
are encoded in the iCS field theory.

While we are interested in the iCS theory as a (3+1)d system, we start by
reviewing the global symmetry and anomaly of a general $\mathrm{U}(1)^{N_{z}}$
CS theory, which can be thought of as compactifying the 3D system in the $z$
direction (so there are $N_{z}$ “layers”). i.e. it reduces to a quasi-2D
theory for small compactification radius or by allowing non-locality along the
compactified direction.

Notice that below we use differential form notations extensively.

###  III.1 Global symmetry and anomaly

We consider the action of a $\mathrm{U}(1)^{N}$ CS theory:

$\mathcal{S}=\int_{M_{3}}\frac{K_{IJ}}{4\pi}a_{I}\wedge da_{J},$ (37)

where $M_{3}$ is a closed three manifold. $\mathbf{K}$ is a symmetric integer
matrix. First we consider the theory with the Chern-Simons term only and
assume that there are no matter fields. This is a valid assumption since all
vortices are gapped, so well below the gap we can just study the pure gauge
theory.

####  III.1.1 Global symmetry

We now enumerate all the unitary symmetries, including 0-form and 1-form, of
the topological action.

First consider the pure topological action Eq. (37). The theory has the
following discrete symmetries:

$a_{I}\rightarrow W_{IJ}a_{J},$ (38)

where $\mathbf{W}$ belongs to GL$(N,\mathbb{Z})$, i.e. an invertible $N\times
N$ integral matrix, such that
$\mathbf{W}^{\mathsf{T}}\mathbf{K}\mathbf{W}=\mathbf{K}$. Physically, when
$\mathbf{K}$ is non-degenerate, $\mathbf{W}$ corresponds to a permutation of
anyon types in the (2+1)d Abelian topological phase.

In addition, a U(1)N gauge theory also has U(1)N 0-form symmetries, the
conservation of magnetic fluxes. The currents of the symmetries are $\star
da_{I}$, where $\star$ is the Hodge star operator. When $\mathbf{K}$ is non-
degenerate, these symmetries are embedded into the 1-form symmetry group.

We now discuss the 1-form symmetry of the model, first assuming that
$\mathbf{K}$ is invertible. A general 1-form transformation takes the
following form:

$a_{I}\rightarrow a_{I}+q_{I}\lambda,$ (39)

where $\lambda$ is a properly normalized flat connection, i.e. $d\lambda=0$,
and the holonomy of $\lambda$ along any 1-cycle is an integer multiple of
$2\pi$. $q_{I}$ is an arbitrary real number at this point. When
$q_{I}\in\mathbb{Z}$, the transformation can be viewed as a (possibly
singular) gauge transformation. Additionally, one can prove that in order for
the partition function of the theory to be invariant, the $q_{I}$’s have to
satisfy the following quantization condition:

$K_{IJ}q_{I}\in\mathbb{Z}.$ (40)

The derivation can be found in Appendix C. If $K$ is invertible, it implies
that $q_{I}$’s must take the values

$q_{I}=(K^{-1})_{IJ}k_{J},k_{J}\in\mathbb{Z},$ (41)

and the 1-form symmetry forms a finite Abelian group, the determinant group of
$\mathbf{K}$. Each 1-form transformation corresponds to an integer vector
$\mathbf{m}$, modulo those that can be written as $\mathbf{K}\mathbf{l}$ for
$\mathbf{l}\in\mathbb{Z}^{N}$, and is in fact one-to-one correspondent to
anyon types. This is well-known, as the 1-form symmetries are generated by
Wilson loops. The 1-form symmetry is generally anomalous if the corresponding
anyon is not a boson [52].

When $\mathbf{K}$ is degenerate, the null space of $\mathbf{K}$ needs to be
considered separately. Suppose that $\mathbf{v}\in\mathbb{Z}^{N}$ is a null
vector, i.e. $\mathbf{K}\mathbf{v}=0$. Without loss of generality we can
assume $\mathbf{v}$ is primitive (i.e. $\mathrm{gcd}(v_{1},v_{2},\dots)=1$).
In that case, the following 1-form transformation

$\mathbf{a}\rightarrow\mathbf{a}+\alpha\mathbf{v}\lambda$ (42)

is always an exact symmetry of the 4d action (without any $2\pi\mathbb{Z}$
shift), where $\alpha$ is an arbitrary real number in $[0,1)$. Therefore each
null vector of $K$ gives rise to a U(1) 1-form symmetry group. Physically, a
null vector corresponds to a mode of the gauge fields without CS term, which
can be dualized to a superfluid, which has a U(1) 1-form symmetry in the
absence of vortices [7]. More generally, if the dimension of the null space is
$r$, then we have $\mathrm{U}(1)^{r}$ 1-form symmetry.

Now we consider what happens when Maxwell terms are included. The
$\mathrm{U}(1)^{N}$ and 1-form symmetries are not affected at all by the
Maxwell terms. The Maxwell terms may break the 0-form symmetry, however. We
then have only the subgroup that preserves the Maxwell terms as well.

We make connections between the symmetries identified in the CS gauge theory,
especially with the K matrix given by Eq. (30), and those found in the coupled
wire lattice model which are discussed in Sec II.1. The magnetic
$\mathrm{U}(1)$ 0-form symmetry of the $I$-th layer is nothing but the planar
$\mathrm{U}(1)_{\varphi}$ symmetry. The U(1) “1-form” symmetry defined in Eq.
(42) can be identified with the conservation of total vorticity Eq. (14).
Finally, as we will see shortly, we can identify the string order parameter
with a “string” monopole operator.

####  III.1.2 ’t Hooft anomaly

The emergent symmetry group defined in the previous section generally has ’t
Hooft anomalies. The discrete part of the 0-form symmetry group has to be
analyzed on a case-by-case basis. The anomaly of the finite 1-form symmetry
group was studied in Ref. [52], and since it is not particularly relevant to
us we do not go into details. So in the following we study the mixed anomaly
between the 0-form $\mathrm{U}(1)^{N}$ and the 1-form symmetry group
$\mathcal{A}\times\mathrm{U}(1)^{r}$.

We consider, without any loss of generality, a U(1) 0-form symmetry with
charge vector $\mathbf{t}$. In other words, the U(1) current is given by

$j=\sum_{I}t_{I}\star\\!da_{I}.$ (43)

We turn on a general background gauge field $A$ for the U(1) symmetry:

$\frac{t_{I}}{2\pi}\int a_{I}\wedge dA.$ (44)

To check whether this coupling can be compatible with the 1-form symmetries,
perform a general 1-form gauge transformation $a_{I}\rightarrow
a_{I}+q_{I}\lambda_{I}$ where now $\lambda_{I}$ is allowed to be non-flat. The
action changes by $-\frac{q_{I}t_{I}}{2\pi}\int A\wedge d\lambda$. To restore
gauge invariance, the additional term can be cancelled by inflow from a 4d
action

$S_{\text{bulk}}=\frac{q_{I}t_{I}}{2\pi}\int_{M_{4}}B\wedge dA,$ (45)

where $B$ is the background 2-form gauge field.

A similar argument works for the $\mathrm{U}(1)^{r}$ part of the 1-form
symmetry group: under the 1-form gauge transformation $a_{I}\rightarrow
a_{I}+\alpha{v}_{I}\lambda$, to maintain gauge invariance we need to have a
bulk theory given by

$S_{\text{bulk}}=\frac{\mathbf{v}\cdot\mathbf{t}}{2\pi}\int_{M_{4}}B\wedge
dA,$ (46)

where now $B$ is a 2-form U(1) gauge field, transforming under the 1-form
gauge transformation as $B\rightarrow B+\alpha d\lambda$.

###  III.2 The filling anomaly

We now analyze the filling anomaly in the layered CS gauge theory. We are
mainly interested in K matrices given in Eq. (30), although we present the
analysis in a form that applies to more general cases. Such theories have a
discrete translation symmetry, generated by

$T_{z}:a_{I}\rightarrow a_{I+1},$ (47)

provided that the K matrix satisfies $K_{IJ}=K_{I+1,J+1}$. This is the formal
way of identifying the index $I$ as labeling layers in the $z$ direction. We
assume that $\mathbf{K}$ is “short-ranged” in the $z$ direction, which means
that there exists an integer $d\geq 0$ such that $K_{IJ}=0$ if $|I-J|>d$. For
now we assume that either $N_{z}\rightarrow\infty$, or periodic boundary
condition is imposed so $a_{I+N_{z}}=a_{I}$. This way, the layered CS theory
describes a highly anisotropic 3D system.

We are going to show that the (3+1)d theory can exist at any filling. Instead
of directly computing the 3D filling anomaly, we take a detour and study the
compactified system, with periodic boundary condition imposed along the $z$
direction. Suppose that the 3D system has a filling factor $\nu$ (i.e. the
average U(1) charge is $\nu$ per unit cell). When viewed as a quasi-2D system,
the filling factor becomes $N_{z}\nu$. Here the dependence on $N_{z}$ reflects
the 3D nature.

Before going to the details, we briefly review the theory of the filling
anomaly. It has been understood now that even though the filling condition
does not correspond to a true, quantized ’t Hooft anomaly (essentially because
the filling factor is a continuous quantity), its impact on the low-energy
physics can be described in the same theoretical framework. We first present a
formal argument, following the approach in [53]. The standard method to detect
the ’t Hooft anomaly of a global symmetry is to couple the system to the
background gauge field of the symmetry, and compute the topological response.
While translation is a spatial symmetry, at low energy it can be described
effectively as an internal symmetry. To this end, we formally introduce
$\mathbb{Z}$ gauge fields $x^{1}$ and $x^{2}$ for the 2D translation symmetry.
We also turn on a U(1) background gauge field $A$. The filling condition can
be captured by the following response

$S_{\text{bulk}}=\nu\int dA\cup x^{1}\cup x^{2}.$ (48)

Before getting to the CS gauge theory, we show how the filling anomaly can be
embedded into the $\mathrm{U}(1)^{[0]}\times\mathrm{U}(1)^{[1]}$ ’t Hooft
anomaly discussed in Sec. III.1. Formally, we start from the bulk action

$\frac{k}{2\pi}\int B\cup dA.$ (49)

To realize the filling anomaly, we demand that the translation background
gauge field is activated through the 2-form background:

$B=2\pi\nu x^{1}\cup x^{2}.$ (50)

Then we see that the action gives a filling factor $k\nu$.

We give a physical interpretation as follows. Eq. (50) implies that when a
loop charged under the U(1)[1] symmetry moves, it picks up a phase factor
given by the flux of the background 2-form gauge field $B$ through the
transverse area swept by the loop. Now recall that $x^{1,2}$ are translational
gauge fields, $x^{1}\cup x^{2}$ can be heuristically interpreted as the area
in the 2D plane. So Eq. (50) basically says that a phase $2\pi\nu\mathcal{A}$,
where $\mathcal{A}$ is the enclosed area, is attached to the loop. The anomaly
action then implies that the loop is identified with magnetic field lines of
the $\mathrm{U}(1)^{[0]}$ gauge field. Therefore, the anomaly action and the
particular background Eq. (50) together can be summarized by the following
intuitive picture: a $2\pi$ U(1)[0] flux picks up $e^{2\pi i\nu}$ phase factor
when it moves around a unit area. This is exactly the Aharonov-Bohm phase
expected from the filling factor.

Such a mechanism is realized by a (2+1)d superfluid. A $2\pi$
$\mathrm{U}(1)^{[0]}$ magnetic flux is trapped in a superfluid vortex. It is
well-known that under boson-vortex duality, vortices see a background magnetic
field, the strength of which is fixed by the filling factor.

We can apply these considerations to the CS theory. The anomaly is given by
Eq. (46). The charge vector is simply $t_{I}=1$, as we have shown in the
duality mapping. We posit that the translation symmetry is realized via the
$\mathrm{U}(1)$ 1-form symmetry corresponding to the “zero mode”
$\mathbf{v}=(1,1,\cdots,1)^{\mathsf{T}}$, by the 2-form background given in
Eq. (50). Together, we find that the coefficient of the topological action is
$\mathbf{v}\cdot\mathbf{t}\nu=N_{z}\nu$. which is the expected filling factor.
In fact, it is clear that as long as the charge vector $\mathbf{t}$ and
$\mathbf{v}$ are translation-invariant (i.e. invariant under, or a certain
multiple), the theory can be defined at arbitrary filling. Physically,
whenever a zero mode exists, after compactification the system becomes a
quasi-2D superfluid where the U(1)[0] symmetry is spontaneously broken.

Notice that so far we have assumed periodic boundary condition in the $z$
direction. Now we turn to open boundary condition. In this case, the $z$
translation symmetry is broken. The filling of the quasi-2D system is still
$N_{z}\nu$. The boundary condition at the top and bottom layers now becomes
crucial. We have considered two kinds of boundary conditions in the coupled
wire model in Sec. II.4.1, which lead to different K matrices after duality
transformation. If the boundary condition preserves all the symmetries, we
found the following K matrix:

$\mathbf{K}_{\text{sym,OBC}}\propto\begin{pmatrix}-1&1&&&0\\\ 1&-2&1&&\\\
&\ddots&\ddots&\ddots&\\\ &&1&-2&1\\\ 0&&&1&-1\end{pmatrix}.$ (51)

With this K matrix, $(1,1,\dots,1)^{\mathsf{T}}$ is still a zero mode, and the
same argument, as used for the periodic boundary condition, applies.

We also studied a symmetry-breaking surface in Sec. II.4.1, which yields the
following K matrix:

$\mathbf{K}_{\text{broken,OBC}}\propto\begin{pmatrix}-2&1&&&0\\\ 1&-2&1&&\\\
&\ddots&\ddots&\ddots&\\\ &&1&-2&1\\\ 0&&&1&-2\end{pmatrix}.$ (52)

This K matrix is non-degenerate, thus representing a fully gapped 2D phase. As
discussed in [43] and [52], such a CS theory has a finite emergent symmetry
group, incompatible with a generic filling. Thus the only way that the theory
can emerge is to to break the symmetries explicitly. This is indeed the case
in the coupled wire model.

###  III.3 Monopole operators

U(1) gauge theories generally admit monopole (or “disorder”) operators, which
are charged under the magnetic U(1) symmetries. For $\mathrm{U}(1)^{N}$ gauge
group, a monopole operator is labeled by the magnetic charge vector
$\mathbf{m}\in\mathbb{Z}^{N}$, i.e. $2\pi m_{I}$ flux for the $a_{I}$ gauge
field. Physically, it is actually an instanton that inserts $2\pi\mathbf{m}$
flux at a point. Due to the Chern-Simons coupling, the flux insertion
necessarily nucleates charges, so we expect that the monopole operator is not
gauge-invariant. In order to build a gauge-invariant monopole operator, one
has to dress the “bare” operator by matter fields of charge
$\mathbf{K}\mathbf{m}$. As matter fields are absent (or very massive) in the
low-energy theory, a gauge-invariant monopole operator can be introduced at
low energy (i.e. well below the mass gap) only when $\mathbf{m}$ is a null
vector of $\mathbf{K}$. For a translation-invariant K matrix, null vectors are
also translation-invariant. They are precisely the string order parameters in
the coupled wire model. Via the well-known Polyakov mechanism, adding monopole
operators to the theory causes confinement of the gauge theory. However, if
the monopole operator has a translation-invariant magnetic charge vector, it
is in fact a highly non-local, string-like object in 3D and cannot appear as a
local term in the Hamiltonian. We thus conclude that such theories are stable
with respect to confinement caused by monopole proliferation.

We now look closer at the monopole operators, which requires adding Maxwell
terms in the action:

$S_{\mathrm{Maxwell}}=\int\tilde{g}_{IJ}f_{I}\wedge\star f_{J}.$ (53)

For this purpose, we diagonalize the K matrix:

$\sum_{J}K_{IJ}e_{J}^{\alpha}=\lambda^{q}e_{I}^{q},$ (54)

where $q$ labels different eigenvectors. Here we normalize $e^{q}$ such that
$\sum_{I}e^{q}_{I}e^{q^{\prime}}_{I}=\delta^{qq^{\prime}}$. Then the K matrix
can be written as $K_{IJ}=\sum_{q}\lambda^{q}e^{q}_{I}e^{q}_{J}$. Define new
gauge fields $b^{q}=e_{I}^{q}a_{I}$ and the inverse transformation reads
$a_{I}=e^{q}_{I}b^{q}$. In terms of the new fields, the CS term becomes

$\sum_{q}\frac{\lambda^{q}}{4\pi}b^{q}\wedge db^{q}.$ (55)

and the Maxwell term becomes

$\tilde{g}_{IJ}e^{q}_{I}e^{q^{\prime}}_{J}f^{q}\wedge\star f^{q^{\prime}},$
(56)

here $f=db$.

For the simplest choice, we can set $\tilde{g}_{IJ}=\tilde{g}\delta_{IJ}$,
then
$\tilde{g}_{IJ}e^{q}_{I}e^{q^{\prime}}_{J}=\tilde{g}\delta^{qq^{\prime}}$. So
the theory becomes

$\mathcal{L}=\frac{\lambda^{q}}{4\pi}b^{q}\wedge
db^{q}-\tilde{g}f^{q}\wedge\star f^{q}.$ (57)

In our layered CS theory, the translation invariance along $z$ guarantees that
different modes do not couple. In general, $\tilde{g}^{a}$ stays finite even
when $N\rightarrow\infty$, unless $\tilde{g}_{IJ}$ is exactly proportional to
$K_{IJ}$.

As we have already mentioned, the gauge-invariant monopole operators
correspond to null space of the $\mathbf{K}$ matrix. Such operators only exist
when the K matrix is degenerate, so it has zero eigenvalues, the number of
which is the same as the dimension of the null space. We focus on one such
mode, and denote by $b^{0}$. Since the CS term vanishes, the $b^{0}$ mode is
just a pure Maxwell theory dual to a superfluid. The charged operator of the
superfluid is the monopole operator.

A $2\pi m_{I}$ instanton for $a_{I}$ corresponds to a $2\pi e_{I}^{q}m_{I}$
instanton for the $b^{q}$ gauge fields. At this point, we pause to discuss the
Dirac quantization condition. For $a^{I}$ gauge fields, we have the standard
quantization condition:

$\int_{M_{2}}\frac{da_{I}}{2\pi}\in\mathbb{Z},$ (58)

where $M_{2}$ is a closed surface. In terms of the new field strength:

$\sum_{\alpha}e^{q}_{I}\int_{M_{2}}\frac{f^{q}}{2\pi}\in\mathbb{Z}.$ (59)

Since we are only interested in gauge-invariant operators, we can assume that
$\int_{M_{2}}f^{q}=0$ for $\lambda^{q}\neq 0$. For simplicity, we assume that
$\mathbf{K}$ has a unique zero mode $\mathbf{v}$. The normalized eigenvector
is $\mathbf{e}^{0}=\frac{\mathbf{v}}{\absolutevalue{\mathbf{v}}}$. Then the
quantization condition is

$\frac{v_{I}}{\absolutevalue{\mathbf{v}}}\int_{M_{2}}\frac{f^{0}}{2\pi}\in\mathbb{Z}.$
(60)

Since gcd$(v_{1},v_{2},\dots)=1$, we find

$\int_{M_{2}}\frac{f^{0}}{2\pi\absolutevalue{\mathbf{v}}}\in\mathbb{Z}.$ (61)

This is the formal way of saying that the fundamental monopole is labeled by
$\mathbf{v}$ (when expressed using the original gauge fields $a_{I}$).

To restore the standard normalization, we need to rescale
$f^{0}\rightarrow\absolutevalue{\mathbf{v}}f^{0}$, so the Lagrangian density
of the Maxwell term of $b^{0}$ now reads

$\mathcal{L}_{0}=-\absolutevalue{\mathbf{v}}^{2}\tilde{g}f^{0}\wedge\star
f^{0}.$ (62)

This theory can be dualized to a superfluid, but now the superfluid density is
proportional to $\frac{1}{\absolutevalue{\mathbf{v}}^{2}\tilde{g}}$. So the
mass generated by instanton proliferation is
$\propto\frac{1}{|\mathbf{v}|^{2}}$ (there is a further fugacity factor, which
should actually decay exponentially with $\absolutevalue{\mathbf{v}}^{2}$).

Now we specialize to the K matrix of the compactified system. In that case, we
have $|\mathbf{v}|=\sqrt{N_{z}}$, thus the superfluid density goes down as
$N_{z}^{-1}$. The monopole operator, which is charged under the magnetic
symmetry of all the layers, is indeed the string superfluid order parameter of
the coupled wire model identified in Sec. II.3. The $N_{z}^{-1}$ decay of the
superfluid density also agrees with the result of the two-point correlation
function of the order parameter.

One notices an apparent discrepancy between the superfluid density which
decays as $N_{z}^{-1}$ and a finite superconducting response. The resolution
is that the charge of the order parameter grows with $N_{z}$ to compensate the
decay of the superfluid density.

###  III.4 Classification of charge excitations

We now discuss the topological classification of charges in the Chern-Simons
theory. Again we consider the case of a finite number $N_{z}$ of layers.
Suppose the null space of the K matrix is spanned by null vectors
$\mathbf{v}_{j},j=1,2,\dots,r$ where $\mathbf{Kv}_{j}=0$. As shown earlier,
each null vector corresponds to a 1-form $\mathrm{U}(1)^{(j)}$ symmetry, and
the corresponding magnetic U(1) 0-form symmetry is spontaneously broken. One
can then show that the total vorticity (with respect to the weak superfluid of
the $\mathrm{U}(1)^{(j)}$ symmetry) is given by
$\mathbf{v}_{j}^{\mathsf{T}}\mathbf{l}$, by e.g. checking the U(1)(j) 1-form
charge of the corresponding Wilson loop.

One can also consider the equivalence classes of non-vortex gauge charges. The
definition of equivalence is essentially the same as that of a topological
Abelian CS theory with a non-degenerate K matrix, which we briefly review now:
a gauge charge $\mathbf{l}$ is a local excitation if and only if it can be
written as $\mathbf{l}=\mathbf{Kl}^{\prime}$ for some integer vector
$\mathbf{l}^{\prime}$. This is true even when $\mathbf{K}$ is degenerate, as
already mentioned in Sec. III.3, and easily follows from the equation of
motion. Two charges are equivalent if their difference is local. From this
definition, it is easy to see that all local excitations must have zero
vorticity, as expected in any superfluid. We thus focus on the non-vortex
charges, the equivalence classes of which form a finite Abelian group, which
is the generalization of the anyon group. To determine the structure of this
group, we can simply perform a $\mathrm{GL}(N_{z},\mathbb{Z})$ transformation
to decouple the null space. In other words, we can find an invertible integer
matrix $\mathbf{W}$, such that

$\mathbf{W}^{\mathsf{T}}\mathbf{KW}=\begin{pmatrix}\mathbf{0}_{r\times r}&\\\
&\tilde{\mathbf{K}}\end{pmatrix},$ (63)

where $\tilde{\mathbf{K}}$ is a non-degenerate matrix. The structure of the
anyon group is determined by the Smith normal form of $\tilde{K}$. The number
of distinct equivalence classes is given by $|\det\tilde{\mathbf{K}}|$. For
the K matrix arising from the coupled wire construction, we find that
$\tilde{\mathbf{K}}$ is basically the same as the K matrix with symmetry-
breaking OBC, but with $N_{z}-1$ layers, and
$|\det\tilde{\mathbf{K}}|=m^{N_{z}-1}\cdot N_{z}$. These results all agree
very well that those of the coupled wire model.

##  IV Summary and discussion

In this work we have introduced a new type of compressible matter: a “weak
superfluid”, which is a 3D anisotropic phase distinct from the conventional
superfluid or Fermi liquid. The weak symmetry breaking is characterized by a
“rod” order parameter, which is supported on a straight line along a fixed
direction, and at the same time there exist no local order parameters. We
argue that the low-energy physics can be captured by a gapless infinite-
component Chern-Simons theory introduced in Ref. [30].

We determine the emergent symmetry of the highly anisotropic “weak
superfluid”, both from the coupled wire model and from the iCS field theory.
Most importantly, the weak symmetry breaking leads to a “cylindrical” U(1)
1-form symmetry, where the symmetry transformations are defined on cylinders
along the $z$ direction (with arbitrary shape in the $xy$ plane), as depicted
in Fig 2. The mixed anomaly between the cylindrical U(1) 1-form symmetry and
the U(1) 0-form symmetry is responsible for the arbitrary filling condition,
i.e. being a compressible phase.

Recently, Ref. [8] identified an interesting conceptual relation between zero-
temperature DC conductivity and “fluxibility”, which is a mixed anomaly
between charge U(1) symmetry and an emergent symmetry. The latter is also
responsible for the arbitrary filling condition. This led to the conjecture
[8] that a translation-invariant quantum system is compressible if and only if
it is fluxible. In our work, the filling condition is indeed realized in the
low-energy theory as a mixed anomaly between the U(1) charge conservation
symmetry and the cylindrical 1-form symmetry, which suggests that the system
is indeed “fluxible”, in an appropriate sense. However, the highly anisotropic
emergent symmetry group results in more complicated zero-temperature DC
transport, quite different from the isotropic theories considered in Ref. [8].

Our microscopic model was given by a coupled wire construction which stitches
together arrays of 1D Luttinger Liquids. We can also understand the phase by
imagining gluing together layers of 2D superfluids. In previous works [31]
[32], other classes of 3D coupled wire models were shown to exhibit fractonic
behavior. Interestingly, for appropriate choice of parameters, these models
can also be made compressible (i.e. invariant under the anomalous
$\mathrm{U}(1)\times\mathrm{U}(1)$ symmetry defined in Eq. (11)). Similar to
the models studied in this work, they become gapless and exhibit weak symmetry
breaking, i.e. there are no local order parameter of any kind. But the non-
local order parameter is now supported on a membrane that spans the whole 2D
section perpendicular to the wires. It is not yet clear whether the physics
can be understood in terms of some variation of the iCS theory. We leave a
detailed investigation of these other models for future work.

Another natural question is whether similar WSB phenomena occur in two-
dimensional compressible systems. While we are not aware of any fundamental
reasons prohibiting such phenomena in 2D, it is unlikely that a 2D coupled
wire model can realize it. The duality transformation employed in this work
can be applied to a general coupled wire model, under reasonable conditions.
The result is then a 2+1D Abelian CS theory, which is either fully gapped when
the K matrix is non-degenerate, or dual to a superfluid when the K matrix is
degenerate, which now has a local order parameter (e.g. the quasi-2D limit of
the model studied in this work after compactification). It will be interesting
to understand whether this is just the limitation of the coupled wire
construction, or there is a more fundamental obstruction in 2D.

The iCS models have been shown to realize new and novel phases of matter. The
gapped case discussed in [30] provides new examples of type-I fractonic order
without any foliation structure. Thus far the iCS theories that have been
analyzed have only involved nearest neighbor couplings. It would be
interesting to study models with a wider range of couplings between the layers
of gauge theories, and explore possible connections with other fracton phases.
We have shown here that the coupled wire construction provides a useful tool
for building microscopic models which realize these phases at low energy.

###### Acknowledgements.

M.C. would like to thank X. Chen for enlightening conversations and
collaboration on a previous project, and D. Else and T. Senthil for useful
correspondence and sharing unpublished results. J.S. would like to thank N.
Read, X. Ma and S. Bryant for helpful discussions. J.S. especially thanks D.
Mross and E. Leviatan for explaining subtleties in their recent work [45]. We
are grateful to John McGreevy for his feedback on the manuscript. M.C.
acknowledges support from NSF under award number DMR-1846109 and the Alfred P.
Sloan foundation. A.D. is supported by the Simons Foundation through the
collaboration on Ultra-Quantum Matter (651438, XC) and by the Institute for
Quantum Information and Matter, an NSF Physics Frontiers Center (PHY-1733907).
_Note added:_ We would like to draw the reader’s attention to an upcoming work
on closely related topics by X.-Q. Ma, H.-T. Lam and X. Chen [54].

##  Appendix A String operators for solitonic excitations

In this section we consider the mobility of the solitons. We focus on the most
elementary excitation, namely a unit soliton of a single $\Theta_{\mathbf{r}}$
term. We work in the mean-field approximation, assuming large $vg$.

One can apply a vertex operator to create solitons, and move them around in
the $yz$ plane. In this model, it is easy to see that $e^{i\theta}$ can hop a
1-soliton along $y$. In other words, $\prod_{y}e^{i\theta_{yz}(x)}$ is a
string operator. Note that strictly speaking, the string creates perfectly
sharp solitons, which are only valid at $g\rightarrow\infty$. At finite $g$,
the soliton is smeared over a length scale $\xi$ and the string operator must
be modified to create the smeared profile. But these modifications do not
change universal features, such as mobility or the braiding statistics.

We also consider how a single soliton excitation can be transported along the
wire. For simplicity we consider the $n_{1}=n_{2}=n$ case, but the same method
works for the more general cases as well. Following the general discussion in
Ref. [31], we define

$\varphi_{yz}^{L}=\varphi_{yz}+m\theta_{yz}+n(\theta_{y,z-1}+\theta_{y,z+1}).$
(64)

We construct a string operator of the following form:

$W_{y}(x_{2},x_{1})=\exp\left(i\int_{x_{1}}^{x_{2}}\sum_{z}w_{z}\partial_{x}\varphi^{L}_{yz}\right).$
(65)

where $w_{z}\in\mathbb{R}$. At the end point $x_{2}$, $W$ creates a soliton of
strength $\sum_{z^{\prime}}K_{zz^{\prime}}w_{z^{\prime}}$ for the
$\Theta_{yz}$ term, where $y$ is fixed, and the opposite one at $x_{1}$. We
thus denote such an excitation at one end of the string by a
$N_{z}$-dimensional integer vector $\mathbf{v}$, the $z$-th entry of which is
the strength of the soliton in the $\Theta_{yz}$ term. Thus finding a string
operator of the form Eq. (65) reduces to solving the equation
$\mathbf{Kw}=\mathbf{v}$.

Thus to move a single soliton of unit strength at $z=0$, we set
$w_{z}=(K^{-1})_{0z}$. When the bulk is fully gapped (i.e.
$2|m|>|n_{1}+n_{2}|$), $w$ as a function of $z$ decays exponentially away from
the location of the excitation (it is strictly localized in $y$), so our
string operator is quasi-localized, and the soliton can move along the wire.
This agrees with the prediction of the iCS field theory [30], that is the
quasiparticles have exponentially localized profiles leading to quasi-
localized braiding statistics. In the coupled wire model, the quasiparticle
excitation itself is strictly localized as a violation of a certain plaquette
term (at least in the mean-field limit), but the string operator is quasi-
localized.

In the gapless case, the construction in Eq. (65) does not work in general as
$K$ is singular, at least when PBC along the $z$ direction is imposed.
However, we can still invert the K matrix in the complement of the null space.
For $m=-2n$, this subspace consists of all vectors $\mathbf{v}$ such that
$\sum_{i}v_{i}\neq 0$. For the vectors with non-zero overlap with the null
space, as we discussed in the main text they should be thought of as vortices
in the weak superfluid. For example, if $\mathbf{v}_{0}=(1,-1,0,\cdots)$, then
we can find the following $\mathbf{w}$ which satisfies
$\mathbf{Kw}=\mathbf{v}_{0}$:

$\mathbf{w}=\frac{1}{m}\Big{(}-\\!\frac{1}{L_{z}},\frac{L_{z}-2}{L_{z}},\frac{L_{z}-1}{L_{z}},\cdots,\frac{1}{L_{z}},0\Big{)}.$
(66)

$\mathbf{w}$ is no longer quasi-localized. The weight is only inversely
proportional to the distance away from the location of the excitation. A
string operator moving this excitation along the wire direction can then be
constructed using the $\mathbf{w}$.

It is shown in Ref. [32] that two excitations $\mathbf{v}$ and
$\mathbf{v}^{\prime}$ (with the same $x$ coordinate) belong to the same
superselection sector, i.e. they can be transformed into each other by acting
with local operators, if and only if
$\mathbf{v}^{\prime}=\mathbf{v}+\mathbf{Ku}$ for some integer vector
$\mathbf{u}$. This agrees with the mathematical definition of superselection
sectors in Abelian CS theories. Let us now enumerate the number of
superselection sectors. We can easily show that, all translations of
$m\mathbf{v}_{0}$ are equivalent, and $mL_{z}\mathbf{v}_{0}$ is local. It is
then straightforward to show that the total number of inequivalent non-vortex
excitations is $m^{L_{z}-1}\cdot L_{z}$. More generally, we can find the Smith
normal form of $\mathbf{K}$, and the number of superselection sectors is the
same as the absolute value of the product of all non-zero entries. Notice that
here we only count the “non-vortex” excitations, as there are infinite types
of vortices labeled by vorticity.

Finally, we consider the mobility of excitations along $z$. Since
$e^{i\theta}$ only moves solitons along $y$, motion along $z$ necessarily
involves $e^{i\varphi}$. It is useful to represent the charge configuration
formally by a Laurent polynomial in two variables $y$ and $z$. Namely, define

$\sum_{i,j\in\mathbb{Z}}q_{ij}y^{i}z^{j}.$ (67)

Here $q_{ij}$ is the charge at coordinate $i,j$ in the $y-z$ plane and it
should be clear that we do not need to worry about their $x$ coordinates. Such
polynomial representation was discussed in Ref. [31].

The operator $e^{i\varphi}$ creates the pattern $(nz^{-1}+m+nz)(1+y)$. Since
$1-y$ can be created by $e^{i\theta}$, for mobility along $z$ we can move all
excitations to the same $y$ and thus just consider $2(nz^{-1}+m+nz)$. Suppose
by applying $e^{ia_{j}\varphi_{j}}$ (together with appropriate $e^{i\theta}$’s
to move all charges to the same $y$), we can create a configuration with two
solitons of opposite charges separated by a distance $l$. In terms of the
polynomial representation, this amounts to finding a polynomial
$f(z)=\sum_{j}a_{j}z^{j}$ such that

$f(z)(nz^{-1}+m+nz)\propto 1+z^{l}.$ (68)

Suppose that $l$ is large so we essentially look for a string operator along
$z$. For $1\ll j\ll l$ we have

$na_{j+1}+ma_{j}+na_{j-1}=0.$ (69)

Since $m^{2}>4n^{2}$, the corresponding characteristic polynomial has two real
roots. As a result, $a_{j}$ grows exponentially with $j$ (in either
directions), which means the charge created at one end costs an exponentially
large amount of energy. For $m=-2n$, let
$f(z)=a_{0}+a_{1}z+\cdots+a_{l}z^{l}$, then
$f(z)(nz^{-1}-2n+nz)=n[a_{0}+(a_{1}-2a_{0})z+(a_{2}-2a_{1}+a_{0})z^{2}+\cdots+a_{l}z^{l+1}]$.
If we require $a_{j+1}-2a_{j}+a_{j-1}=0$ for all $1\leq j\leq l-1$, then we
$a_{j}=a_{0}+j(a_{1}-a_{0})$. So the only way to have bounded $|a_{j}|$ is to
set $a_{1}=a_{0}$, and the excitation created at one end is precisely a
vortex-anti-vortex dipole of strength $n$.

##  Appendix B Details of the duality mapping

In this section we provide more detail for the mapping described in section
II.4.

###  B.1 Boson-vortex duality

The wire model is described by an array of Luttinger liquids with conjugate
variables $(\varphi,\theta)$ and Lagrangian

$\begin{split}\mathcal{L}[\varphi,\theta]&=\sum_{\mathbf{r}}\frac{i}{\pi}\partial_{x}\theta_{\mathbf{r}}\partial_{\tau}\varphi_{\mathbf{r}}+\frac{\tilde{v}}{2\pi}\left(\partial_{x}\varphi_{{\mathbf{r}}}\right)^{2}+\frac{u}{2\pi}\left(\partial_{x}\theta_{{\mathbf{r}}}\right)^{2}+\frac{v}{8\pi}\left(\partial_{x}\Delta_{y}\varphi_{{\mathbf{r}}}\right)^{2}-g\cos(2\bar{\theta}_{\mathbf{r}})\end{split}$
(70)

where

$\begin{split}2\bar{\theta}_{{\mathbf{r}}}&=\Delta_{y}\varphi_{\mathbf{r}}+\Lambda\theta_{\mathbf{r}}\\\
&=\left(\varphi_{{\mathbf{r}}+\hat{y}}-\varphi_{{\mathbf{r}}}\right)+\left(m\theta_{{\mathbf{r}}}+m\theta_{{\mathbf{r}}+\hat{y}}+n_{1}\theta_{{\mathbf{r}}-\hat{z}}+n_{2}\theta_{{\mathbf{r}}+\hat{z}}+n_{2}\theta_{{\mathbf{r}}+\hat{y}-\hat{z}}+n_{1}\theta_{{\mathbf{r}}+\hat{y}+\hat{z}}\right).\end{split}$
(71)

Here we define
$\Delta_{y}X_{\mathbf{r}}=X_{{\mathbf{r}}+\hat{y}}-X_{\mathbf{r}}$ and
${\mathbf{r}}=(y,z)$ is a wire index. The definition of $\Lambda$ can be
easily inferred from Eq. (71). Note also that
$\frac{v}{8\pi}\left(\partial_{x}\Delta_{y}\varphi_{{\mathbf{r}}}\right)^{2}$
has been added to the standard Luttinger Liquid kinetic term. This is done
because it is convenient for the duality mapping to the vortex theory but its
presence does not affect the qualitative physics [44].

In what follows we use the two equivalent forms of labelling for the wires,
$O_{{\mathbf{r}}+\hat{a}}\equiv O_{y+a_{y},z+a_{z}}$. Anticipating the layered
structure of the gauge theory, we treat $z$ as the “layer” index and $y$ is to
be coarse grained to a continuous spatial coordinate. With this motivation in
mind, we define the following pair of conjugate variables:

$\begin{array}[]{l}\tilde{\varphi}_{yz}=-\sum_{y^{\prime}}\operatorname{sgn}\left(y^{\prime}-y-\frac{1}{2}\right)\theta_{y^{\prime}z}\\\
\\\
\tilde{\theta}_{yz}=\frac{1}{2}\left(\varphi_{y+1,z}-\varphi_{y,z}\right).\end{array}$
(72)

Here $\tilde{\varphi}_{yz}$ creates a $2\pi$ vortex in the $\varphi$ field in
layer $z$, in between wires $y$ and $y+1$. The operator
$\partial_{x}\tilde{\theta}_{yz}$ is the “charge” operator for this vortex. We
can re-express the Lagrangian in Eq. (70) in terms of these new fields but the
result will be highly non-local in the $y$ direction. We can restore locality
via a Hubbard-Stratonovich transformation at the expense of introducing new
degrees of freedom $a_{0}^{(z)}(x,y)$ and $a_{1}^{(z)}(x,y)$ in each layer
indexed by $(z)$.

We now rewrite the kinetic part of the Lagrangian using the dual variables
$(\tilde{\varphi},\tilde{\theta})$ and the Hubbard-Stratonovich fields
$(a_{0}^{(z)},a_{1}^{(z)})$:

$\begin{array}[]{l}\sum_{\mathbf{r}}\frac{i}{\pi}\partial_{x}\theta_{\mathbf{r}}\partial_{\tau}\varphi_{\mathbf{r}}=\sum_{\mathbf{r}}\frac{i}{\pi}\partial_{x}\tilde{\theta}_{\mathbf{r}}\partial_{\tau}\tilde{\varphi}_{\mathbf{r}}\\\
\\\
\sum_{\mathbf{r}}\frac{v}{8\pi}\left(\partial_{x}\Delta_{y}\varphi_{{\mathbf{r}}}\right)^{2}=\sum_{\mathbf{r}}\frac{v}{2\pi}\left(\partial_{x}\tilde{\theta}_{{\mathbf{r}}}\right)^{2}\\\
\\\
\sum_{\mathbf{r}}\frac{\tilde{v}}{2\pi}\left(\partial_{x}\varphi_{{\mathbf{r}}}\right)^{2}=\sum_{\mathbf{r}}\frac{\tilde{v}}{2\pi}\left(\Delta_{y}^{-1}\partial_{x}\tilde{\theta}_{\mathbf{r}}\right)^{2}\to\sum_{\mathbf{r}}\left[-\frac{i}{\pi}\partial_{x}\tilde{\theta}_{\mathbf{r}}a_{0}^{(z)}+\frac{(\Delta_{y}a_{0}^{(z)})^{2}}{8\pi\tilde{v}}\right]\\\
\\\
\sum_{\mathbf{r}}\frac{u}{2\pi}\left(\partial_{x}\theta_{{\mathbf{r}}}\right)^{2}=\sum_{\mathbf{r}}\frac{u}{8\pi}\left(\partial_{x}\Delta_{y}\tilde{\varphi}_{{\mathbf{r}}}\right)^{2}\rightarrow\sum_{\mathbf{r}}\frac{u}{8\pi}\left[\left(\partial_{x}\tilde{\varphi}_{{\mathbf{r}}}-a_{1}^{(z)}\right)^{2}+\left(\Delta_{y}a_{1}^{(z)}\right)^{2}\right]+\sum_{z,y,y^{\prime}}\frac{u}{8\pi}V_{y,y^{\prime}}\partial_{x}\left(\Delta_{y}\tilde{\varphi}_{yz}\right)\partial_{x}\left(\Delta_{y}\tilde{\varphi}_{y^{\prime}z}\right)\end{array}$
(73)

where
$V=\Delta_{y}^{\mathsf{T}}(1+\Delta_{y}^{\mathsf{T}}\Delta_{y})^{-1}\Delta_{y}$,
which decays exponentially in the difference in wire index $y$. Up to this
point everything is exact. Now we interpret $a_{\mu}^{(z)}$ as a (2+1)d gauge
field, living in layer $z$, in the $a_{2}^{(z)}=0$ gauge. The upshot of all
this is that we get a theory of layers of vortices minimally coupled to a
layered gauge theory. Using Eq. (72), we re-express the plaquette term in
terms of the vortex variables

$\begin{split}2\bar{\theta}_{yz}&=2\tilde{\theta}_{yz}+\frac{1}{2}\Delta_{y}\left[n_{1}(\tilde{\varphi}_{y,z+1}+\tilde{\varphi}_{y+1,z-1})+n_{2}(\tilde{\varphi}_{y,z-1}+\tilde{\varphi}_{y+1,z+1})+m(\tilde{\varphi}_{yz}+\tilde{\varphi}_{y+1,z})\right]\\\
&=2\tilde{\theta}_{yz}+\frac{1}{2}\Lambda\cdot\Delta_{y}\tilde{\varphi}_{yz}~{}~{}.\end{split}$
(74)

In what follows, when coupled to matter fields $(\bar{\varphi},\bar{\theta})$,
we suppress the $(z)$ layer index of the gauge fields and use the
${\mathbf{r}}$ index in order to make the expressions more succinct. The
various indices can be parsed as follows
$a_{\mu,{\mathbf{r}}}=a_{\mu,yz}=a_{\mu}^{(z)}(x,y)$. At this stage the
Lagrangian takes the form

$\mathcal{L}[\tilde{\varphi},\tilde{\theta},a_{\mu}]=\sum_{\mathbf{r}}\frac{i}{\pi}\partial_{x}\tilde{\theta}_{\mathbf{r}}\left(\partial_{\tau}\tilde{\varphi}_{\mathbf{r}}-a_{0,{\mathbf{r}}}^{(z)}\right)+\frac{u}{2\pi}\left(\partial_{x}\tilde{\varphi}_{{\mathbf{r}}}-a_{1,{\mathbf{r}}}^{(z)}\right)^{2}+\frac{v}{2\pi}\left(\partial_{x}\tilde{\theta}_{\mathbf{r}}\right)^{2}-g\cos(2\tilde{\theta}_{\mathbf{r}}+\frac{1}{2}\Lambda\cdot\Delta_{y}\tilde{\varphi}_{\mathbf{r}})+\mathcal{L}_{\text{Maxwell}}^{(z)}.$
(75)

Here, in the $a_{2}^{(z)}=0$ gauge,
$\mathcal{L}_{\text{Maxwell}}^{(z)}=\frac{u}{2\pi}(\Delta_{y}a_{1}^{(z)})^{2}+\frac{1}{2\pi\tilde{v}}(\Delta_{y}a_{0}^{(z)})^{2}$.
In order to simplify the sine-Gordon term we can introduce another set of
conjugate variables:

$\bar{\varphi}_{{\mathbf{r}}}=\tilde{\varphi}_{{\mathbf{r}}}~{}~{}~{}\text{and}~{}~{}~{}\bar{\theta}_{\mathbf{r}}=\tilde{\theta}_{\mathbf{r}}+\frac{1}{4}\Lambda\cdot\Delta_{y}\bar{\varphi}_{\mathbf{r}}.$
(76)

Expressing the Lagrangian in this final basis and also expanding the cosine
term yields

$\begin{split}\mathcal{L}[\bar{\varphi},\bar{\theta},a_{\mu}]=&\sum_{\mathbf{r}}\frac{i}{\pi}\partial_{x}\bar{\theta}_{\mathbf{r}}(\partial_{\tau}\bar{\varphi}_{\mathbf{r}}-a_{0,{\mathbf{r}}})+\frac{i}{4\pi}\left(\Lambda\cdot\Delta_{y}\partial_{x}\bar{\varphi}\right)_{\mathbf{r}}a_{0,{\mathbf{r}}}+\frac{u}{2\pi}\left(\partial_{x}\bar{\varphi}_{{\mathbf{r}}}-a_{1,{\mathbf{r}}}\right)^{2}\\\
&~{}~{}~{}+\frac{v}{2\pi}\left(\partial_{x}\bar{\theta}_{\mathbf{r}}\right)^{2}-\frac{v}{4\pi}\left(\Lambda\cdot\Delta_{y}\partial_{x}\bar{\varphi}\right)_{\mathbf{r}}\partial_{x}\bar{\theta}_{\mathbf{r}}+\frac{v}{32\pi}\left(\Lambda\cdot\Delta_{y}\partial_{x}\bar{\varphi}\right)_{\mathbf{r}}^{2}\\\
&~{}~{}~{}+g\bar{\theta}_{\mathbf{r}}^{2}+\mathcal{L}_{\text{Maxwell}}.\end{split}$
(77)

###  B.2 Integrating out matter fields

In order to get a pure gauge theory the plaquette degrees of freedom
$\bar{\varphi},\bar{\theta}$ must be integrated out. The first step in the
derivation is to collect all the terms involving $\bar{\varphi}$, up to
$\mathcal{O}(\bar{\varphi}^{2})$:

$\begin{split}&\mathcal{O}(\bar{\varphi}^{2}):~{}\partial_{x}\bar{\varphi}_{{\mathbf{r}}^{\prime}}\left[\frac{u}{2\pi}+\frac{v}{32\pi}(\Lambda\cdot\Delta_{y})^{\mathsf{T}}(\Lambda\cdot\Delta_{y})\right]_{{\mathbf{r}}^{\prime}{\mathbf{r}}}\partial_{x}\bar{\varphi}_{\mathbf{r}}=\partial_{x}\bar{\varphi}M\partial_{x}\bar{\varphi}\\\
&\mathcal{O}(\bar{\varphi}):~{}\underbrace{\left[\frac{-i}{4\pi}\Lambda\cdot\Delta_{y}a_{0}-\frac{u}{\pi}a_{1}+\frac{i}{\pi}\partial_{\tau}\bar{\theta}+\frac{v}{4\pi}\Lambda\cdot\Delta_{y}\partial_{x}\bar{\theta}\right]_{\mathbf{r}}}_{\Gamma_{\mathbf{r}}}\partial_{x}\bar{\varphi}_{\mathbf{r}}\end{split}$
(78)

We denote the cross term by
$\Gamma_{\mathbf{r}}[\bar{\theta},a]\partial_{x}\bar{\varphi}_{\mathbf{r}}.$
Note here we have used a discrete “integration by parts” in parts involving
$\Lambda\cdot\Delta_{y}$ and dropped $z$ boundary terms. This step requires
care when considering open boundary conditions in $z$. Integrating out
$\bar{\varphi}$ results in an expression of the form:

$-\frac{1}{4}\Gamma_{{\mathbf{r}}^{\prime}}\left(M^{-1}\right)_{{\mathbf{r}}^{\prime}{\mathbf{r}}}\Gamma_{\mathbf{r}}.$
(79)

We are interested in the low energy theory so we can express $M^{-1}$ as a
derivative expansion in powers of $\Lambda\cdot\Delta_{y}$. To accomplish this
note that for $M=1+D$

$\begin{split}M^{-1}&=\left(1+D\right)^{-1}=\sum_{j=0}^{\infty}(-1)^{j}D^{j}~{}.\end{split}$
(80)

In the present case
$M^{-1}\approx\frac{2\pi}{u}-\frac{v\pi}{8u^{2}}(\Lambda\cdot\Delta_{y})^{\mathsf{T}}(\Lambda\cdot\Delta_{y})$.
We keep the 2nd order term because the $a_{1}^{2}$ term in $\Gamma^{2}$
cancels and the next lowest order term is $(\Lambda\cdot\Delta_{y}a_{1})^{2}$.
Including the left over terms which do not involve $\bar{\varphi}$ and only
keeping the lowest order derivative terms in the expansion of $\Gamma^{2}$
gives

$\begin{split}\mathcal{L}[\bar{\theta},a]&=\sum_{\mathbf{r}}\frac{-i}{4\pi}(\Lambda\cdot\Delta_{y}a_{0})_{\mathbf{r}}a_{1,{\mathbf{r}}}+\frac{v}{8\pi}(\Lambda\cdot\Delta_{y}a_{1})_{\mathbf{r}}^{2}+\frac{1}{32\pi
u}(\Lambda\cdot\Delta_{y}a_{0})_{\mathbf{r}}^{2}+\frac{i}{\pi}\bar{\theta}_{\mathbf{r}}\left(\partial_{x}a_{0}-\partial_{\tau}a_{1}\right)_{\mathbf{r}}+g\bar{\theta}_{\mathbf{r}}^{2}\\\
&~{}~{}~{}~{}+\mathcal{L}_{\text{Maxwell}}+\text{Higher order
terms}\end{split}$ (81)

Integrating out $\bar{\theta}$ is now relatively simple and the result is

$\begin{split}\mathcal{L}[a]=&\sum_{\mathbf{r}}-\frac{i}{4\pi}\underbrace{(\Lambda\cdot\Delta_{y}a_{1})_{\mathbf{r}}a_{0,{\mathbf{r}}}}_{\text{CS
term}}+\frac{1}{4\pi^{2}g}(\partial_{x}a_{0}-\partial_{\tau}a_{1})_{\mathbf{r}}^{2}\\\
&~{}~{}~{}~{}~{}~{}~{}+\frac{1}{2\pi\tilde{v}}(\Delta_{y}a_{0})_{\mathbf{r}}^{2}+\frac{1}{32\pi
u}(\Lambda\cdot\Delta_{y}a_{0})_{\mathbf{r}}^{2}+\frac{u}{2\pi}(\Delta_{y}a_{1})_{\mathbf{r}}^{2}+\frac{v}{8\pi}(\Lambda\cdot\Delta_{y}a_{1})_{\mathbf{r}}^{2}.\end{split}$
(82)

If the $y$ direction is coarse grained and the layer $(z)$ index restored then

$\begin{split}\left(\Lambda\cdot\Delta_{y}a\right)_{z}&\approx(n_{1}+n_{2})\Delta_{y}a^{(z+1)}+2m\Delta_{y}a^{(z)}+(n_{1}+n_{2})\Delta_{y}a^{(z-1)}\end{split}$
(83)

So the first term in Eq. (82), corresponding to a CS term, becomes

$\mathcal{L}_{\text{CS}}[a]=\frac{i}{4\pi}\left[(n_{1}+n_{2})a_{1}^{(z+1)}\partial_{y}a_{0}^{(z)}+2ma_{1}^{(z)}\partial_{y}a_{0}^{(z)}+(n_{1}+n_{2})a_{1}^{(z-1)}\partial_{y}a_{0}^{(z)}\right]~{}~{}~{}.$
(84)

When the gauge constraint, $a_{2}=0$, is relaxed
$a^{(z)}_{1}\partial_{y}a_{0}^{(z^{\prime})}\to\epsilon^{\mu\nu\lambda}a^{(z)}_{\mu}\partial_{\nu}a_{\lambda}^{(z^{\prime})}$,
so the CS part of the Lagrangian has the more familiar form

$\mathcal{L}_{\mathrm{CS}}[a]\equiv\frac{iK_{zz^{\prime}}}{4\pi}a^{(z)}\wedge
da^{(z^{\prime})}~{}~{}.$ (85)

We can similarly investigate parts of Eq. (82) associated with the Maxwell
terms of the gauge theory. For example, coarse graining $y$ and restoring
gauge invariance in the $4^{\text{th}}$ term corresponds to

$(\Lambda\cdot\Delta_{y}a_{0})^{2}\to(\Lambda\cdot[\Delta_{y}a_{0}-\partial_{\tau}a_{2}])^{2}=\left[(n_{1}+n_{2})E_{2}^{(z+1)}+2mE_{2}^{(z)}+(n_{1}+n_{2})E_{2}^{(z-1)}\right]^{2}~{}~{}.$
(86)

We see that the typical $E_{y}^{2}$ term present in conventional Maxwell
theory is replaced by something which couples the same component ($y$ in this
case) of the electric field in different layers. Indeed, a similar term is
generated for $B$ but not for $E_{x}$, this is not surprising given the
anisotropic nature of the underlying microscopic model. The full result for
the Maxwell term is given by

$\mathcal{L}_{\text{Maxwell}}=\sum_{zz^{\prime}}E_{1}^{(z)}\frac{\delta_{zz^{\prime}}}{4\pi^{2}g}E_{1}^{(z^{\prime})}+E_{2}^{(z)}\left[\frac{\delta_{zz^{\prime}}}{2\pi\tilde{v}}+\frac{1}{32\pi
u}(K^{\mathsf{T}}K)_{zz^{\prime}}\right]E_{2}^{(z^{\prime})}+B^{(z)}\left[\frac{u\delta_{zz^{\prime}}}{2\pi}+\frac{v}{8\pi}(K^{\mathsf{T}}K)_{zz^{\prime}}\right]B^{(z^{\prime})}~{}~{}.$
(87)

The upshot of this discussion is that the effective IR Lagrangian is given by:

$\mathcal{L}[a]\equiv\frac{iK_{zz^{\prime}}}{4\pi}a^{(z)}\wedge
da^{(z^{\prime})}+\mathcal{L}_{\text{Maxwell}}~{}~{}.$ (88)

Note that for the compressible state, $m=-(n_{1}+n_{2})$,
$K^{\mathsf{T}}K\to\partial_{z}^{4}$ under coarse-graining in $z$ and so the
inter-layer portion of the Maxwell term would appear to be irrelevant in this
context.

###  B.3 Coupling to vortices

So far we have assumed that $\varphi$ and $\theta$ vary smoothly, excluding
the “solitonic” excitations of the cosine pinning term. These more singular
configurations are the analog of vortices in the usual boson-vortex duality,
which should be minimally coupled to the dynamical gauge fields. we now extend
the duality map to include these excitations.

To this end, we replace the cosine potential
$g\cos(2\bar{\theta}_{\mathbf{r}})$ by

$\frac{-g}{2}(2\bar{\theta}_{\mathbf{r}}-2\pi n_{\mathbf{r}})^{2}$ (89)

Here $n_{\mathbf{r}}$ parametrizes the locations of the solitons where
$2\bar{\theta}$ jumps by integer multiples of $2\pi$.

From the discussion of the microscopic model we can view $n$ as corresponding
to some configuration of vortices created by an operator of the form
$e^{i\sum\int n(x,y,z)\partial_{x}\tilde{\varphi}}$. It is natural to wonder
how the introduction of such a field manifests in the dual gauge theory. This
can be settled using the mapping of section II.4. Indeed, amending Eq 81 gives

$\begin{split}\mathcal{L}[\bar{\theta},a]&=\sum_{\mathbf{r}}\frac{-i}{4\pi}(\Lambda\cdot\Delta_{y}a_{0})_{\mathbf{r}}a_{1,{\mathbf{r}}}+\frac{v}{8\pi}(\Lambda\cdot\Delta_{y}a_{1})_{\mathbf{r}}^{2}+\frac{1}{32\pi
u}(\Lambda\cdot\Delta_{y}a_{0})_{\mathbf{r}}^{2}\\\
&~{}~{}~{}~{}+\frac{i}{\pi}\bar{\theta}(\partial_{x}a_{0}-\partial_{\tau}a_{1}+2\pi
ign)+g\bar{\theta}_{\mathbf{r}}^{2}+\frac{g}{4}n_{\mathbf{r}}^{2}\\\
&~{}~{}~{}~{}+\mathcal{L}_{\text{Maxwell}}+\text{Higher order
terms}.\end{split}$ (90)

Integrating out $\bar{\theta}$ gives Eq. (82) plus a coupling term between the
gauge fields and $n$

$\mathcal{L}[a,n]=\mathcal{L}[a]+i(\partial_{x}a_{0}^{(z)}-\partial_{\tau}a_{1}^{(z)})n\rightarrow\mathcal{L}[a]-(a_{0}i\partial_{x}n-a_{1}i\partial_{\tau}n)\\\
.$ (91)

This tells us that we can interpret $n$ as a $x$-dipole of the gauge charge.
Namely, $\partial_{x}n=\rho$ is the charge density. One can expect that if we
treat the dynamics of solitons carefully (e.g. solitons moving between wires),
one would obtain the full matter-gauge coupling $j\cdot a$. Therefore, as
expected, vortices become gauge charges.

##  Appendix C Quantization condition for 1-form symmetry transformation

In order to derive the quantization conditions on $q_{I}$, it is more
convenient to use the 4d definition of the Chern-Simons term. Suppose that
$M_{3}$ is the boundary of a 4d manifold $M_{4}$, i.e. $\partial M_{4}=M_{3}$.
The gauge field is also extended to $M_{4}$. Then the CS action can be defined
as

$S=\frac{K_{IJ}}{4\pi}\int_{M_{4}}F_{I}\wedge F_{J}.$ (92)

Now consider 1-form transformations parametrized by $\lambda_{I}$. We assume
that they are extended to the bulk as well, but not necessarily flat. The 4d
action changes:

$\begin{split}\delta
S&=\frac{K_{IJ}}{4\pi}\int_{M_{4}}(q_{I}d\lambda_{I}\wedge
F_{J}+q_{J}F_{I}\wedge d\lambda_{J}+q_{I}q_{J}d\lambda_{I}\wedge
d\lambda_{J})\\\ &=\frac{K_{IJ}}{4\pi}\int_{M_{4}}(q_{I}d\lambda_{I}\wedge
F_{J}+q_{J}F_{I}\wedge
d\lambda_{J})+\frac{K_{IJ}}{4\pi}q_{I}q_{J}\int_{M_{3}}\lambda_{I}\wedge
d\lambda_{J}\\\ &=\frac{K_{IJ}}{4\pi}\int_{M_{4}}(q_{I}d\lambda_{I}\wedge
F_{J}+q_{J}F_{I}\wedge d\lambda_{J})\\\ &=2\pi
K_{IJ}q_{I}\int_{M_{3}}\frac{\lambda_{I}}{2\pi}\wedge\frac{F_{J}}{2\pi}\end{split}$
(93)

Since $\int_{M_{3}}\frac{\lambda_{I}}{2\pi}\wedge\frac{F_{J}}{2\pi}$ is an
integer (intersection number), for the partition function to remain the same,
we must have $e^{i\delta S}=1$, or $K_{IJ}q_{I}$ an integer.

##  Appendix D Phenomenologies

###  D.1 Order parameter

The stretched superfluid order parameter is given by the string operator
$e^{i\Phi_{y}}$ where $\Phi_{y}(x)=\sum_{z}\varphi_{yz}$. We want to compute
the two point function: $\langle
e^{i\Phi_{y}(x)}e^{-i\Phi_{y^{\prime}}(x^{\prime})}\rangle$. To tackle this we
approximate the sine-Gordon term by the quadratic term $\Theta^{2}$ in the
strong coupling limit. The resulting theory is free so we can shift our
attention to calculating the
$\langle\Phi_{y}(x)\Phi_{y^{\prime}}(x^{\prime})\rangle$, which is equal to

$\begin{split}\sum_{z,z^{\prime}}\langle\varphi_{yz}(x)\varphi_{y^{\prime}z^{\prime}}(x^{\prime})\rangle&=\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{q}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{q^{\prime}}\langle\varphi_{q}\varphi_{q^{\prime}}\rangle\sum_{z,z^{\prime}}e^{iq\cdot
r}e^{iq^{\prime}\cdot
r^{\prime}}=\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}G_{\varphi}(q)e^{ik_{\perp}\cdot(r-r^{\prime})}\sum_{z,z^{\prime}=1}^{N_{z}}e^{ik_{z}(z-z^{\prime})}.\end{split}$
(94)

Here $k_{\perp}=(\omega,k_{x},k_{y})$ is the momentum perpendicular to $z$,
$q=(k_{\perp},k_{z})$ and $G(q,q^{\prime})=\delta_{q,q^{\prime}}G(q)$, since
the theory is free and translation invariant. The sum over $z,z^{\prime}$ can
be easily evaluated

$\sum_{zz^{\prime}}e^{ik_{z}(z-z^{\prime})}=N_{z}\sum_{z}e^{ik_{z}z}=N_{z}\delta(k_{z}).$
(95)

Thus

$\begin{split}\langle\Phi_{y}(x)\Phi_{y^{\prime}}(x^{\prime})\rangle&=N_{z}\int
dk_{\perp}\,e^{ik_{\perp}\cdot{\mathbf{r}}}G_{\varphi}(k_{z}=0).\end{split}$
(96)

To work out the two-point function $G_{\varphi}$ we re-express the action as

$S=\sum_{q}\begin{pmatrix}\varphi(q)&\theta(q)\end{pmatrix}G^{-1}(q)\begin{pmatrix}\varphi(-q)\\\
\theta(-q)\end{pmatrix}~{}.$ (97)

Recall the mean-field Lagrangian Eq. (25) is given by:

$S=\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{{\mathbf{r}},x,\tau}dxd\tau\frac{i}{\pi}\partial_{x}\theta_{\mathbf{r}}\partial_{\tau}\varphi_{\mathbf{r}}+\left[\tilde{v}(\partial_{x}\varphi_{\mathbf{r}})^{2}+u(\partial_{x}\theta_{\mathbf{r}})^{2}+g\Theta_{\mathbf{r}}^{2}\right]~{}~{}.$
(98)

Using this, we get

$G^{-1}(q)=\begin{pmatrix}\tilde{v}k_{x}^{2}+gf_{\varphi}(q)&\frac{ik_{x}\omega}{\pi}+gf(q)\\\
&\\\ \frac{ik_{x}\omega}{\pi}+gf(-q)&uk_{x}^{2}+gf_{\theta}(q)\end{pmatrix}$
(99)

where

$\begin{array}[]{l}f_{\varphi}=2-2\cos k_{y}\\\ \\\
f_{\theta}=4\left[m\cos\frac{k_{y}}{2}+n_{1}\cos(\frac{k_{y}}{2}-k_{z})+n_{2}\cos(\frac{k_{y}}{2}+k_{z})\right]^{2}\\\
\\\ f=2i\big{[}n_{1}\sin k_{z}+n_{2}\sin(k_{y}+k_{z})+m\sin k_{y}+n_{1}\sin\
(k_{y}-k_{z})-n_{2}\sin k_{z}\big{]}~{}~{}~{}.\\\ \end{array}$ (100)

We are interested in realizations with the chiral U(1) present which occurs
when $n_{1}+n_{2}=-m$. In this case $f_{\varphi}f_{\theta}-|f|^{2}=0$ and so
we can express $G_{\varphi}$ as follows

$G_{\varphi}=\frac{uk_{x}^{2}+gf_{\theta}}{k_{x}^{2}\left[\tilde{v}uk_{x}^{2}+\frac{\omega^{2}}{\pi^{2}}+g\left(uf_{\varphi}+\tilde{v}f_{\theta}\right)\right]}.$
(101)

Setting $k_{z}=0$, we have $f_{\theta}=0$, and $G_{\varphi}$ is simplified to

$G_{\varphi}=\frac{u}{\frac{\omega^{2}}{\pi^{2}}+u\tilde{v}k_{x}^{2}+2gu(1-\cos
k_{y})}.$ (102)

This form makes sense as $k_{z}=0$ means $\theta$ fields are independent of
$z$, so the wire coupling term decouples into
$\cos(\varphi_{yz}-\varphi_{y+1,z})$. For long wavelength limit, we can expand
$1-\cos k_{y}\approx\frac{k_{y}^{2}}{2}$, and rescale
$k_{x}\rightarrow\frac{k_{x}}{\sqrt{\tilde{v}}},k_{y}\rightarrow\frac{k_{y}}{g}$,

$\langle\Phi_{y}(x)\Phi_{0}(0)\rangle\sim\frac{N_{z}}{\sqrt{gx^{2}+\tilde{v}y^{2}}},$
(103)

which can be understood as a superfluid with anisotropic dispersion. Thus the
two-point function for our order parameter has the following scaling:

$\langle e^{i\Phi_{y}(x)}e^{-i\Phi_{0}(0)}\rangle\sim
e^{-C\frac{N_{z}}{\sqrt{gx^{2}+\tilde{v}y^{2}}}}.$ (104)

Here $C$ is a constant.

###  D.2 Compressibility

The presence of both global U(1) symmetries means the model can be defined at
arbitrary filling. This, in turn, suggests compressibility. Here we add more
field theoretic justification for this feature. To compute compressibility we
add a uniform chemical potential term to the Hamiltonian

$\begin{split}H-\mu
N&=H-\mu\sum_{\mathbf{r}}\int\partial_{x}\theta_{\mathbf{r}}\\\
&=\int\frac{1}{2\pi}[\tilde{v}(\partial_{x}\varphi)^{2}+u(\partial_{x}\theta)^{2}]-\mu\partial_{x}\theta-g\cos(2\bar{\theta})\end{split}$
(105)

where we have suppressed $\sum_{\mathbf{r}}$ in the last line. The term
$\mu\partial_{x}\theta$ can be absorbed into $(\partial_{x}\theta)^{2}$ by
completing the square:
$\frac{u}{2\pi}(\partial_{x}\theta)^{2}-\mu\partial_{x}\theta=\frac{u}{2\pi}(\partial_{x}\theta-\frac{\pi}{u}\mu)^{2}-\frac{\pi}{2u}\mu^{2}$.
Next we can shift $\partial_{x}\theta-\frac{\pi}{u}\mu\to\partial_{x}\theta$
by redefining $\theta\to\theta+\frac{\pi}{u}\mu x$ but this is also a symmetry
of the sine-Gordon term. In particular under this transformation
$2\bar{\theta}\to 2\bar{\theta}$. Including the uniform chemical potential the
action is given by

$S[\mu]=S-\frac{\pi V}{2u}\mu^{2}.$ (106)

The compressibility is $\kappa\sim\frac{\partial^{2}F}{\partial\mu^{2}}$
which, for this action, is finite. Note that the above argument rests on the
fact that the shift $\theta\to\theta+\frac{\pi}{u}\mu x$ is a symmetry of
$2\bar{\theta}$. This is a consequence of the presence of the $U(1)_{\theta}$
subsystem symmetry. So we are in some way factoring in the stretched order
parameter in this calculation.

###  D.3 Electromagnetic response

In order to work out the transport properties of the model we can consider
coupling the gauge theory to a 3+1D background field $A$. The natural current
is the monopole current in each layer
$j_{\mu}^{I}\sim\epsilon^{\mu\nu\lambda}\partial_{\nu}a^{I}_{\lambda}$. Indeed
the appropriate minimal coupling can be derived from the microscopic model by
replacing $\partial_{\mu}\varphi\to\partial_{\mu}\varphi+A_{\mu}$ and carrying
out the mapping of section II.4. For simplicity consider the case with an
intra-layer Maxwell term and a common charge $t$ in each layer as in section
III.3,

$\mathcal{L}=\frac{K^{IJ}}{4\pi}a_{I}\wedge da_{J}-\tilde{g}f_{I}\wedge\star
f_{I}+\frac{it}{2\pi}A_{I}\wedge da_{I}$ (107)

where here $A_{I}=A(x,y;z=I)$. Note the $A_{z}(x,y,z)$ does not couple to any
current, which reflects the fact that the transport of charge along layers is
highly suppressed. This can be seen from the form of the coupling term in the
microscopic model. Recall the sine-Gordon term
$\Theta_{yz}=\varphi_{yz}-\varphi_{y+1,z}+(\theta\text{ terms })$.
Heuristically this condenses a process in which a charged particle created by
$e^{i\varphi}$ is hopped along the y direction with some additional
complicated back-scattering of vortices between different $z$ layers, but no
charges tunnel between layers. As was done in Sec. III.3, we can diagonalize
Eq. (107).

$\mathcal{L}=\sum_{q}\frac{\lambda^{q}}{4\pi}b^{q}\wedge
db^{q}-\tilde{g}f^{q}\wedge\star f^{q}+\frac{it}{2\pi}A^{q}\wedge db^{q}.$
(108)

In what follows we use $q$ subscripts and superscripts interchangeably, (e.g.
$O_{q}\equiv O^{q}$) since there is no subtlety with raising and lowering
operators for the layer index.

It is helpful to integrate out the gauge fields $b$ in order to get a theory
purely in terms of the background field $A$. Note that
$\mathcal{L}=\sum_{q}\mathcal{L}_{q}$ where, defining
$m_{q}=\frac{\lambda_{q}}{8\pi q}$ and working in the
$\partial_{\mu}b^{q}_{\mu}=0$ gauge,

$\begin{split}\mathcal{L}_{q}&=-\tilde{g}f^{\mu\nu,q}f_{\mu\nu}^{q}+2m_{q}\epsilon^{mu\nu\gamma}b_{\mu}^{q}\partial_{\nu}b^{q}_{\gamma}+b_{\mu}^{q}J^{\mu,q}=2\tilde{g}~{}b^{q}_{\mu}\left(\partial^{2}\eta^{\mu\nu}-m_{q}\epsilon^{\gamma\mu\nu}\partial_{\gamma}\right)b^{q}_{\nu}+b_{\mu}^{q}J^{\mu,q}~{}~{}~{}.\end{split}$
(109)

One can check that the propagator for this theory is given by

$\begin{split}D^{\mu\nu}_{qq^{\prime}}(x,y)=\frac{-\delta_{qq^{\prime}}}{2\tilde{g}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{k^{2}-m_{q}^{2}}\left[\eta^{\mu\nu}-m_{q}^{2}\frac{k^{\mu}k^{\nu}}{k^{4}}+i\frac{m_{q}}{k^{2}}\epsilon^{\mu\nu\gamma}k_{\gamma}\right]e^{ik\cdot(x-y)}~{}~{}.\end{split}$
(110)

Integrating out $b_{q}$ results in the effective action
$S[J]=-\frac{1}{2}\iint J(x)D(x,y)J(y)$. Here the current is given by
$J^{\mu,q}=\frac{it}{2\pi}\epsilon^{\mu\nu\gamma}\partial_{\nu}A_{\gamma}^{q}$.
Working in momentum space the full effective action is given by

$\begin{split}S[A]&=\frac{t^{2}}{16\pi^{2}\tilde{g}}\epsilon^{\mu\alpha\gamma}\epsilon^{\nu\beta\xi}\sum_{q}\int\frac{d^{3}k}{(2\pi)^{3}}k_{\alpha}A^{q}_{\gamma}(k)D^{qq}_{\mu\nu}(k)k_{\beta}A_{\xi}(-k)\\\
&=\frac{t^{2}}{16\pi^{2}\tilde{g}}\sum_{q}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{k^{2}(A^{q})^{2}-\left(k\cdot
A^{q}\right)^{2}-im_{q}\epsilon^{\mu\nu\gamma}A^{q}_{\mu}(k)k_{\nu}A^{q}_{\gamma}(-k)}{k^{2}-m_{q}^{2}}~{}~{}.\end{split}$
(111)

Restricting attention to the $q=0$ sector, where $m_{q}=0$, the Lagrangian has
the form

$\mathcal{L}_{q=0}=\frac{t^{2}}{16\pi^{2}\tilde{g}}A_{\mu}^{q=0}(k)\left(\eta^{\mu\nu}-\frac{k^{\mu}k^{\nu}}{k^{2}}\right)A^{q=0}_{\nu}(-k)~{}~{}.$
(112)

From this we can surmise that the application of a long wavelength ($q\to
0,k\to 0$ ) vector potential induces a Meissner effect, a hallmark of
superconductivity [55].

###  D.4 Interaction between vortices

Making the $(z)$ layer index $I$ of the background field manifest, we have
seen that vortices are minimally coupled to the gauge field:

$\begin{split}\mathcal{L}[a,n]&=\mathcal{L}[a]+\sum_{I}a^{I}_{\mu}j^{\mu,I},\end{split}$
(113)

where $j^{\mu,I}$ is the vortex current in the $I$-th layer and we have used
the layered Maxwell-CS theory of section III.3:

$\mathcal{L}[a]=\frac{K_{IJ}}{4\pi}a_{I}\wedge
da_{J}+\tilde{g}\sum_{z}f_{I}\wedge\star f_{I}.$ (114)

As before we can diagonalize the K-matrix. For the compressible model i.e. for
$m=-(n_{1}+n_{2})$, the normalized eigenvectors
$e^{q}_{I}=\frac{1}{\sqrt{N}}e^{iqI}$ have eigenvalues $\lambda_{q}=2-2\cos
q$. In this basis $\mathcal{L}=\sum_{q}\mathcal{L}_{q}$ with

$\begin{split}\mathcal{L}_{q}=\frac{\lambda_{q}}{4\pi}b^{q}\wedge
db^{q}+\tilde{g}f^{q}\wedge\star f^{q}+b^{q}\wedge\star j^{q}.\end{split}$
(115)

Integrating out the gauge fields, we find the effective action of the currents
$j_{\mu}$ is [56]

$E=\frac{1}{2T}\sum_{qq^{\prime}}\iint
d^{3}xd^{3}y\,j^{\mu}_{q}(x)D_{\mu\nu}^{qq^{\prime}}(x,y)j^{\nu}_{q^{\prime}}(y).$
(116)

where the propagator $D^{\mu\nu}_{qq^{\prime}}$ is given by Eq. (110). Here
the measure is given by
$d^{3}x=dx^{0}dx^{1}dx^{2}=dtdxdy=dtd^{2}{\mathbf{r}}$. Notice that here
${\mathbf{r}}$ denotes the position in the $xy$ plane, different from our
convention in the main text.

Consider a vortex-antivortex pair at $({\mathbf{r}}_{a},z_{a})$ and
$({\mathbf{r}}_{b},z_{b})$ respectively. This corresponds to
$n(x,y;z)=2\pi\delta_{zz_{a}}\delta_{yy_{a}}\Theta(x-x_{a})-2\pi\delta_{zz_{b}}\delta_{yy_{b}}\Theta(x-x_{b})$
which can be expressed as the following current

$j^{\mu}_{(z)}=i\eta^{\mu}_{0}[\delta({\mathbf{r}}-{\mathbf{r}}_{a})\delta_{zz_{a}}-\delta({\mathbf{r}}-{\mathbf{r}}_{b})\delta_{zz_{b}}]=\sum_{q}i\eta^{\mu}_{0}[\delta({\mathbf{r}}-{\mathbf{r}}_{a})e^{-iqz_{a}}-\delta({\mathbf{r}}-{\mathbf{r}}_{b})e^{-iqz_{b}}]\frac{e^{iqz}}{\sqrt{N_{z}}}=\sum_{q}j^{\mu}_{q}\frac{e^{iqz}}{\sqrt{N_{z}}}$
(117)

which minimally couples to the gauge field
$j^{\mu}_{(z)}a^{(z)}_{\mu}=j^{\mu,I}a_{\mu}^{I}$.

Plugging this into Eq 116 and ignoring the self energy terms gives the
following

$\begin{split}E&=\frac{1}{2T}\sum_{qq^{\prime}}\iint
d^{3}xd^{3}y\frac{e^{iqz_{a}}}{\sqrt{N_{z}}}\frac{e^{iq^{\prime}z_{b}}}{\sqrt{N_{z}}}D_{qq^{\prime}}^{00}\delta^{2}(x-a)\delta^{2}(y-b)\\\
&=\frac{-1}{4T\tilde{g}}\sum_{q}\iint
dx^{0}dy^{0}\frac{e^{iq(z_{a}-z_{b})}}{N_{z}}\int\frac{d^{3}p}{(2\pi)^{3}}~{}e^{-ix_{0}p_{0}}e^{-iy_{0}p_{0}}\frac{e^{i{\mathbf{p}}\cdot({\mathbf{r}}_{a}-{\mathbf{r}}_{b})}}{p^{2}-\lambda_{q}^{2}/\tilde{g}^{2}}\left(\eta^{00}-\lambda_{q}^{2}\frac{p^{0}p^{0}}{\tilde{g}^{2}p^{4}}\right)\\\
&=\frac{1}{4TN_{z}\tilde{g}}\sum_{q}\int
dx^{0}e^{iq(z_{a}-z_{b})}\int\frac{dp^{1}dp^{2}}{(2\pi)^{2}}~{}\frac{e^{i{\mathbf{p}}\cdot({\mathbf{r}}_{a}-{\mathbf{r}}_{b})}}{{\mathbf{p}}\cdot{\mathbf{p}}+\lambda_{q}^{2}/\tilde{g}^{2}}\\\
&=\frac{1}{2N_{z}\tilde{g}}\sum_{q}e^{iq(z_{a}-z_{b})}\int\frac{dp^{1}dp^{2}}{(2\pi)^{2}}~{}\frac{e^{i{\mathbf{p}}\cdot({\mathbf{r}}_{a}-{\mathbf{r}}_{b})}}{{\mathbf{p}}\cdot{\mathbf{p}}+\lambda_{q}^{2}/\tilde{g}^{2}}\end{split}$
(118)

where in going from line 2 to line 3 we first used $\int
dx^{0}e^{-ix_{0}p_{0}}=\delta(p_{0})$ and then carried out the integral $\int
dp^{0}$. In going from line 3 to line 4 we used $\int dy^{0}=T$. In the large
$N_{z}$ limit we can replace $\frac{1}{N_{z}}\sum_{q}\to\int dq$, so

$\begin{split}E[r_{a},r_{b}]&\sim\int
d^{2}{\mathbf{p}}dq\frac{e^{i{\mathbf{p}}\cdot({\mathbf{r}}_{a}-{\mathbf{r}}_{b})}e^{iq(z_{a}-z_{b})}}{{\mathbf{p}}^{2}+\frac{1}{\tilde{g}^{2}}(2-2\cos
q)^{2}}\\\ &\sim\int
d^{2}{\mathbf{p}}dq\frac{e^{i{\mathbf{p}}\cdot({\mathbf{r}}_{a}-{\mathbf{r}}_{b})}e^{iq(z_{a}-z_{b})}}{{\mathbf{p}}^{2}+\frac{1}{\tilde{g}^{2}}q^{4}}~{}~{}.\end{split}$
(119)

To evaluate this we identify $q^{2}/\tilde{g}$ with the mass in a 2D free
theory. Placing one of the vortices the origin we have

$\begin{split}E[{\mathbf{r}},z]&\sim\int
dqe^{iqz}K_{0}(q^{2}|{\mathbf{r}}|/\tilde{g})=\frac{\pi^{2}}{\sqrt{8}}\frac{\tilde{g}z}{|{\mathbf{r}}|}\left[I_{-1/4}\left(\frac{z^{2}\tilde{g}}{8|{\mathbf{r}}|}\right)^{2}-I_{1/4}\left(\frac{z^{2}\tilde{g}}{8|{\mathbf{r}}|}\right)^{2}\right].\end{split}$
(120)

Note that in the strong coupling limit $\tilde{g}$ is small.

When $z^{2}/|{\mathbf{r}}|\gg\frac{1}{\tilde{g}}$ then we must use the
$I_{\nu}(x\to\infty)$ asymptotic expansion. Note that
$K_{\nu}=\frac{\pi}{2}\frac{I_{-\nu}-I_{\nu}}{\sin\pi\nu}$, so
$I^{2}_{-\nu}-I^{2}_{\nu}=\frac{2\sin\pi\nu}{\pi}K_{\nu}\left(I_{-\nu}+I_{\nu}\right)$.
Combining the asymptotic expansion for $K_{1/4}$ and $I_{\pm 1/4}$ gives

$E\approx\frac{4\pi}{|z|}+\mathcal{O}\left(\frac{|{\mathbf{r}}|^{2}}{\tilde{g}^{2}z^{5}}\right).$
(121)

In the other limit where $z^{2}/|{\mathbf{r}}|$ is comparable or much less
than 1 we use the asymptotic form of the limit $I_{\nu}(x\to 0)$. In this case
$I^{2}_{-1/4}(u)-I^{2}_{1/4}(u)=\frac{\sqrt{2}}{\Gamma(3/4)^{2}}\frac{1}{\sqrt{u}}-\frac{\sqrt{u}}{\Gamma(5/4)^{2}\sqrt{2}}+\ldots$
and so

$\begin{split}E\approx\sqrt{\frac{\tilde{g}}{|{\mathbf{r}}|}}+\mathcal{O}\left(z^{2}\frac{\tilde{g}^{3/2}}{|{\mathbf{r}}|^{3/2}}\right).\end{split}$
(122)

Having worked out the basic example of the vortex-antivortex pair we can
compute the interaction energy of more complicated configurations of
vorticees. Consider the case of a pair of vortex dipoles; one at the origin
oriented in the $\hat{z}$ direction and one at position
$({\mathbf{r}}_{1},z_{1})$ oriented in the negative $\hat{z}$ direction. The
current is given by

$j^{\mu}_{(z)}=i\eta^{\mu}_{0}[\delta({\mathbf{r}}-{\mathbf{r}}_{1})(\delta_{z,z_{1}}-\delta_{z,z_{1}+1})-\delta({\mathbf{r}})(\delta_{z,0}-\delta_{z,1})].$
(123)

In the limit where $z_{1}^{2}/|{\mathbf{r}}_{1}|\gg\frac{1}{\tilde{g}}$ the
interaction energy is given by

$E\approx\frac{8\pi a_{0}^{2}}{|z_{1}|^{3}},$ (124)

where $a_{0}$ is the distance between neighboring layers. In the limit where
$z_{1}^{2}/|{\mathbf{r}}_{1}|$ is comparable or much less than 1, which
includes the case where the dipoles reside in the same $xy$-plane, the
interaction energy scales as

$E\propto-\left(\frac{\tilde{g}}{|{\mathbf{r}}_{1}|}\right)^{3/2}~{}~{}~{}.$
(125)

## References

  * Ardonne  _et al._ [2004] E. Ardonne, P. Fendley, and E. Fradkin, Topological order and conformal quantum critical points, Ann. Phys. 310, 493–551 (2004).
  * Sur and Yang [2019] S. Sur and K. Yang, Metallic state in bosonic systems with continuously degenerate dispersion minima, Phys. Rev. B 100, 024519 (2019).
  * Lake  _et al._ [2021] E. Lake, T. Senthil, and A. Vishwanath, Bose-luttinger liquids (2021), arXiv:2101.02197 [cond-mat.str-el] .
  * Else  _et al._ [2021] D. V. Else, R. Thorngren, and T. Senthil, Non-fermi liquids as ersatz fermi liquids: General constraints on compressible metals, Phys. Rev. X 11, 021005 (2021).
  * Cheng  _et al._ [2016] M. Cheng, M. Zaletel, M. Barkeshli, A. Vishwanath, and P. Bonderson, Translational symmetry and microscopic constraints on symmetry-enriched topological phases: A view from the surface, Phys. Rev. X 6, 041068 (2016).
  * Cho  _et al._ [2017] G. Y. Cho, C.-T. Hsieh, and S. Ryu, Anomaly manifestation of lieb-schultz-mattis theorem and topological phases, Phys. Rev. B 96, 195105 (2017).
  * Delacrétaz  _et al._ [2020] L. Delacrétaz, D. Hofman, and G. Mathys, Superfluids as higher-form anomalies, SciPost Phys. 8, 10.21468/scipostphys.8.3.047 (2020).
  * Else and Senthil [2021] D. Else and T. Senthil, Critical drag as a mechanism for resistivity (2021), arXiv:2106.15623 [cond-mat.str-el] .
  * Note [1] Similarly, one can consider a stack of 2D Fermi liquids. However, such a state has an open Fermi surface consisting of two curves wrapping around the Brillouin zone.
  * Chamon [2005] C. Chamon, Quantum glassiness in strongly correlated clean systems: An example of topological overprotection, Phys. Rev. Lett. 94, 040402 (2005).
  * Bravyi  _et al._ [2011] S. Bravyi, B. Leemhuis, and B. M. Terhal, Topological order in an exactly solvable 3D spin model, Ann. Phys. 326, 839 (2011), arXiv:1006.4871 [quant-ph] .
  * Haah [2011] J. Haah, Local stabilizer codes in three dimensions without string logical operators, Phys. Rev. A 83, 042330 (2011).
  * Vijay  _et al._ [2015] S. Vijay, J. Haah, and L. Fu, A new kind of topological quantum order: A dimensional hierarchy of quasiparticles built from stationary excitations, Phys. Rev. B 92, 235136 (2015).
  * Vijay  _et al._ [2016] S. Vijay, J. Haah, and L. Fu, Fracton topological order, generalized lattice gauge theory, and duality, Phys. Rev. B 94, 235157 (2016).
  * Williamson [2016] D. J. Williamson, Fractal symmetries: Ungauging the cubic code, Phys. Rev. B 94, 155128 (2016).
  * Nandkishore and Hermele [2019] R. M. Nandkishore and M. Hermele, Fractons, Ann. Rev. Condensed Matter Phys. 10, 295 (2019), arXiv:1803.11196 [cond-mat.str-el] .
  * Pretko  _et al._ [2020] M. Pretko, X. Chen, and Y. You, Fracton Phases of Matter, Int. J. Mod. Phys. A 35, 2030003 (2020), arXiv:2001.01722 [cond-mat.str-el] .
  * Xu [2006] C. Xu, Gapless bosonic excitation without symmetry breaking: An algebraic spin liquid with soft gravitons, Phys. Rev. B 74, 224433 (2006).
  * Rasmussen  _et al._ [2016] A. Rasmussen, Y.-Z. You, and C. Xu, Stable Gapless Bose Liquid Phases without any Symmetry, ArXiv e-prints (2016), arXiv:1601.08235 [cond-mat.str-el] .
  * Pretko [2017a] M. Pretko, Subdimensional particle structure of higher rank $u(1)$ spin liquids, Phys. Rev. B 95, 115139 (2017a).
  * Pretko [2017b] M. Pretko, Generalized electromagnetism of subdimensional particles: A spin liquid story, Phys. Rev. B 96, 035119 (2017b).
  * Pretko [2018] M. Pretko, The fracton gauge principle, Phys. Rev. B 98, 115134 (2018), arXiv:1807.11479 .
  * Seiberg and Shao [2021a] N. Seiberg and S.-H. Shao, Exotic symmetries, duality, and fractons in 2+1-dimensional quantum field theory, SciPost Phys. 10, 10.21468/scipostphys.10.2.027 (2021a).
  * Seiberg and Shao [2020] N. Seiberg and S.-H. Shao, Exotic $U(1)$ Symmetries, Duality, and Fractons in 3+1-Dimensional Quantum Field Theory, SciPost Phys. 9, 046 (2020), arXiv:2004.00015 [cond-mat.str-el] .
  * Seiberg and Shao [2021b] N. Seiberg and S.-H. Shao, Exotic $\mathbb{Z}_{n}$ symmetries, duality, and fractons in 3+1-dimensional quantum field theory, SciPost Phys. 10, 10.21468/scipostphys.10.1.003 (2021b).
  * Slagle and Kim [2017] K. Slagle and Y. B. Kim, Quantum field theory of x-cube fracton topological order and robust degeneracy from geometry, Phys. Rev. B 96, 195139 (2017).
  * Slagle [2021a] K. Slagle, Foliated quantum field theory of fracton order, Phys. Rev. Lett. 126, 101603 (2021a).
  * Shirley  _et al._ [2019] W. Shirley, K. Slagle, and X. Chen, Foliated fracton order from gauging subsystem symmetries, SciPost Physics 6, 10.21468/scipostphys.6.4.041 (2019).
  * You  _et al._ [2018] Y. You, T. Devakul, F. J. Burnell, and S. L. Sondhi, Subsystem symmetry protected topological order, Phys. Rev. B 98 (2018), arXiv:1803.02369 .
  * Ma  _et al._ [2020] X. Ma, W. Shirley, M. Cheng, M. Levin, J. McGreevy, and X. Chen, Fractonic order in infinite-component chern-simons gauge theories (2020), arXiv:2010.08917 [cond-mat.str-el] .
  * Sullivan  _et al._ [2021a] J. Sullivan, A. Dua, and M. Cheng, Fractonic topological phases from coupled wires, Phys. Rev. Research 3, 023123 (2021a).
  * Sullivan  _et al._ [2021b] J. Sullivan, T. Iadecola, and D. J. Williamson, Planar p-string condensation: Chiral fracton phases from fractional quantum hall layers and beyond, Phys. Rev. B 103, 205301 (2021b).
  * You  _et al._ [2020] Y. You, T. Devakul, S. L. Sondhi, and F. J. Burnell, Fractonic chern-simons and bf theories, Phys. Rev. Research 2, 023249 (2020).
  * Williamson  _et al._ [2019] D. J. Williamson, Z. Bi, and M. Cheng, Fractonic matter in symmetry-enriched $u(1)$ gauge theory, Phys. Rev. B 100, 125150 (2019).
  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Ann. Phys. 321, 2–111 (2006).
  * Wang and Levin [2013] C. Wang and M. Levin, Weak symmetry breaking in two-dimensional topological insulators, Phys. Rev. B 88, 245136 (2013).
  * Rao and Sodemann [2021] P. Rao and I. Sodemann, Theory of weak symmetry breaking of translations in z2 topologically ordered states and its relation to topological superconductivity from an exact lattice z2 charge-flux attachment, Phys. Rev. Research 3, 10.1103/physrevresearch.3.023120 (2021).
  * Kane  _et al._ [2002] C. L. Kane, R. Mukhopadhyay, and T. C. Lubensky, Fractional quantum hall effect in an array of quantum wires, Phys. Rev. Lett. 88, 036401 (2002).
  * Teo and Kane [2014] J. C. Teo and C. Kane, From luttinger liquid to non-abelian quantum hall states, Physical Review B 89, 085101 (2014).
  * Fuji and Furusaki [2019] Y. Fuji and A. Furusaki, From coupled wires to coupled layers: Model with three-dimensional fractional excitations, Physical Review B 99, 241107 (2019).
  * Haldane [1995] F. D. M. Haldane, Stability of chiral luttinger liquids and abelian quantum hall states, Phys. Rev. Lett. 74, 2090 (1995).
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Two soluble models of an antiferromagnetic chain, Annals of Physics 16, 407 (1961).
  * Gaiotto  _et al._ [2015] D. Gaiotto, A. Kapustin, N. Seiberg, and B. Willett, Generalized global symmetries, Journal of High Energy Physics 2015, 10.1007/jhep02(2015)172 (2015).
  * Mross  _et al._ [2017] D. F. Mross, J. Alicea, and O. I. Motrunich, Symmetry and duality in bosonization of two-dimensional dirac fermions, Phys. Rev. X 7, 041016 (2017).
  * Leviatan and Mross [2020] E. Leviatan and D. F. Mross, Unification of parton and coupled-wire approaches to quantum magnetism in two dimensions, Phys. Rev. Research 2, 043437 (2020).
  * Imamura  _et al._ [2019] Y. Imamura, K. Totsuka, and T. H. Hansson, From coupled-wire construction of quantum hall states to wave functions and hydrodynamics, Phys. Rev. B 100, 125148 (2019).
  * Slagle [2021b] K. Slagle, Foliated quantum field theory of fracton order, Phys. Rev. Lett. 126, 101603 (2021b).
  * Aasen  _et al._ [2020] D. Aasen, D. Bulmash, A. Prem, K. Slagle, and D. J. Williamson, Topological defect networks for fractons of all types, Phys. Rev. Research 2, 043165 (2020).
  * Wang [2020] J. Wang, Non-liquid cellular states (2020), arXiv:2002.12932 .
  * Wen [2020] X.-G. Wen, Systematic construction of gapped nonliquid states, Phys. Rev. Research 2, 033300 (2020).
  * Note [2] At first glance the spectrum of the layered gauge theory arising from Eq. (87) would seem to depend only on $n_{1}+n_{2}$ and $m$, which does not agree with the spectrum derived for the wire model. This discrepancy stems from the fact the $y$ direction has been coarse grained in the layered picture. Restoring higher order derivative terms in $y$ in the gauge theory gives a spectrum matching Eq. (17).
  * Hsin  _et al._ [2019] P.-S. Hsin, H. T. Lam, and N. Seiberg, Comments on one-form global symmetries and their gauging in 3d and 4d, SciPost Physics 6, 10.21468/scipostphys.6.3.039 (2019).
  * Song  _et al._ [2021] X.-Y. Song, Y.-C. He, A. Vishwanath, and C. Wang, Electric polarization as a nonquantized topological response and boundary luttinger theorem, Phys. Rev. Research 3, 023011 (2021).
  * [54] X. Ma, H. T. Lam, and X. Chen, to appear.
  * Scalapino  _et al._ [1993] D. J. Scalapino, S. R. White, and S. Zhang, Insulator, metal, or superconductor: The criteria, Phys. Rev. B 47, 7995 (1993).
  * Zee [2003] A. Zee,  _Quantum field theory in a nutshell_ (2003).