††institutetext: University of Texas, Austin, Physics Department, Austin TX
78712, USA

# Spontaneously Broken Subsystem Symmetries

Jacques Distler  Andreas Karch  Amir Raz
distler@golem.ph.utexas.edu,karcha<EMAIL_ADDRESS>
###### Abstract

We investigate the spontaneous breaking of subsystem symmetries directly in
the context of continuum field theories by calculating the correlation
function of charged operators. Our methods confirm the lack of spontaneous
symmetry breaking in some of the existing continuum field theories with
subsystem symmetries, as had previously been established based on a careful
analysis of the spectrum. We present some novel continuum field theory
constructions that do exhibit spontaneous symmetry breaking whenever allowed
by general principles. These interesting patterns of symmetry breaking occur
despite the fact that all the theories we study are non-interacting.

††preprint: UTTG-21-2021

##  1 Introduction

Fracton lattice models paramekanti2002ring ; Chamon:2004lew ; haah2011local ;
Nandkishore_2019 ; Pretko_2020  have forced a rethinking of some of the
longest held believes about quantum field theories. Standard lore would posit
that, at very long distances, any quantum field theory should either flow to a
gapped system, potentially described by a non-trivial topological field theory
(TFT), or a gapless system described by a conformal field theory (CFT).
Fractonic lattice models do not fit into either of these frameworks; their
limited mobility excitations as well as their system size dependent
groundstate degeneracy is not describably by any standard TFT or CFT.

A genuine continuum field theory framework capable of describing fractonic
physics has been laid out in Seiberg_2020 ; Seiberg:2020bhn ; Seiberg:2020wsg
. The major new ingredient in these field theories are subsystem symmetries:
global symmetries that do not act on the entire system but only on degrees of
freedom localized along certain sub-manifolds. Such subsystem symmetries are
often present in fractonic lattice models; for the exotic field theories
constructed in Seiberg_2020 ; Seiberg:2020bhn ; Seiberg:2020wsg  they become a
definining feature. While some of these new field theories don’t even exhibit
genuine fractonic excitations, they all have some kind of subsystem symmetry.

At face value these exotic field theories look somewhat conventional. The
simplest examples are free field theories; in fact for many examples the
symmetries prevent any relevant or marginal interactions. On the surface the
only new features are unusual kinetic terms. For example, in the simplest 2+1
dimensional model of Seiberg:2020bhn  based on a single real scalar, the
standard $(\partial_{x}\phi)^{2}+(\partial_{y}\phi)^{2}$ kinetic term gets
replaced with $(\partial_{x}\partial_{y}\phi)^{2}$. This is the lowest
dimension kinetic term consistent with a subsystem shift symmetry demanding
the action be invariant under shifts of $\phi$ that can depend on either $x$
or $y$, $\delta\phi=f_{x}(x)+f_{y}(y)$. This kinetic term, however, leads to
dramatic consequences. Field configurations with arbitrary large
discontinuities in the $x$ direction have low energy, and are hence
unsuppressed in the path integral, as long as they are smooth in $y$, and vice
versa. The presence of these discontinuous field configurations gives rise to
a surprising interplay between UV and IR physics which was originally observed
in Seiberg:2020bhn ; Seiberg:2020wsg , and later expanded on in
you2021fractonic ; you2021fractonic2 ; gorantla2021lowenergy . Most notably,
many quantities in the continuum theory remain sensitive to the lattice
spacing of the underlying discrete model, that is to UV scale physics.

One place where this UV/IR mixing becomes important is when trying to
understand whether the subsystem symmetries are spontaneously broken. Already
in the original examples of Seiberg_2020 ; Seiberg:2020bhn ; Seiberg:2020wsg
the issue of symmetry breaking turned out to be quite subtle. The theory has a
large number of classical zero energy modes, carrying a large momentum in one
direction but no momentum in the other. One could identify these modes with
Goldstone modes of a putatively broken subsystem symmetry. However, quantum
mechanically these modes get pushed to energies of order the UV cutoff. The
symmetry is in fact restored. Closely related is the fact that the lightest
states charged under the subsystem symmetry that could detect any breaking get
pushed to UV scale energies as well.

A very good probe of symmetry breaking are correlation functions of charged
operators. One thing we do in this work is to calculate these correlation
functions directly in the continuum theories of Seiberg_2020 ; Seiberg:2020bhn
; Seiberg:2020wsg  and verify that the subsystem symmetries do, in fact,
remain unbroken as suggested by the spectrum. This same conclusion has also
been recently obtained in gorantla2021lowenergy , where the same correlation
functions were obtained by carefully taking the limit of the discretized
theory.

Given these results, one may wonder whether and when spontaneous breaking of
subsystem symmetries in a continuum field theory can happen. As in ordinary
quantum field theories, where the Mermin-Wagner-Coleman theorem forbids
breaking of a continuous symmetry in 1+1 dimensions and below, spontaneous
symmetry breaking should not always be allowed. A theorem due to Batista and
Nussinov Batista:2004sc  shows that the criteria of whether a subsystem
symmetry acting on a $n$ dimensional submanifold can be broken or not are the
same as for an ordinary symmetry for a field theory living in an $n$
dimensional spacetime. For sufficiently high $n$ there is nothing that forbids
breaking and so one may wonder how to realize this phenomenon in exotic field
theories.

One option to analyze the spontaneous breaking of subsystem symmetries is to
work directly with the lattice models. Indeed, one of the early studies of
fractons showed that the spontaneous breaking of a discrete subsystem symmetry
accounts for the large ground state degeneracy in the X-cube model Vijay_2016
. More recently progress has been made in studying the phase transition
between the broken and unbroken phase of discrete subsystem symmetries in
lattice models you2021fractonic2 ; PhysRevB.98.035112 ; You_2020 ; Qi_2021 ;
lake2021subdimensional ; Bulmash_2018 . The breaking of the continuous
subsystem symmetry of Seiberg:2020bhn ; Seiberg:2020wsg  can also be achieved
by taking a different limit of the lattice model, namely an infinte volume
limit at finite lattice spacing gorantla2021lowenergy . While this limit does
give rise to a broken phase, it does not seem to be amendable to continuum
field theory language.

In contrast to these lattice constructions, in this work we present some novel
continuum field theories that do exhibit subsystem symmetry breaking.
Interestingly enough, all these theories are free theories, just like the
original work of Seiberg_2020 ; Seiberg:2020bhn ; Seiberg:2020wsg , and so our
results are exact. Additionally, we find that whenever the theorem of
Batista:2004sc  does allow symmetry breaking, our models indeed do lead to
spontaneously broken subsystem symmetries.

The outline of this work is as follows. In section 2 we will review why we
chose correlations of charged operators to detect symmetry breaking and show
that our methods reproduce the known results for the theories of Seiberg_2020
; Seiberg:2020bhn ; Seiberg:2020wsg , consistent with the recent paper
gorantla2021lowenergy . In section 3 we introduce the classical XY-plaquette
model as one example of a field theory that does in fact lead to a
spontaneously broken symmetry in dimensions where this is not prohibited by
the theorem of Batista:2004sc . In section 4 we generalize our model to a
large class of closely related theories that all follow the same pattern. We
end with a discussion and future directions in section 5.

##  2 Diagnosing the spontaneous breaking of subsystem symmetries

Let us consider a field theory living in $d$ spacetime dimensions with a
subsystem symmetry acting on $n$ dimensional submanifolds. According to the
modified Mermin-Wagner-Coleman theorem Batista:2004sc , the rules for whether
spontaneous breaking of these subsystem symmetries can or can not occur are
equivalent to the standard considerations for a quantum field theory living in
$n$ dimensions. In particular, spontaneous breaking of a continuous global
subsystem symmetry can only happen if $n\geq 2$. However this theorem in and
of itself does little to predict if a specific subsystem symmetry is
spontaneously broken, or even how to diagnose if such a breaking occurs.

The classical approach to diagnosing a spontaneously broken symmetry is to
define a local order parameter that is charged under the symmetry, and then
calculate if this order parameter acquires an expectation value in the vacuum
state. This computation is delicate as the expectation value of a charged
operator vanishes for any finite size system which obeys the symmetry. Hence
the correct procedure is to calculate the expectation value of the order
parameter in the presence of a symmetry breaking background field, which is
then taken to zero. 111For a modern mathematical description of spontaneous
symmetry breaking see strocchi2005symmetry .

The canonical example of this procedure is to consider a complex field in $d$
dimensions, $\phi$, and a Lagrangian with a Mexican hat potential

$S=\int
d^{d}x\left[\partial_{\mu}\phi^{\dagger}\partial^{\mu}\phi+g\left(\phi^{\dagger}\phi-v^{2}\right)^{2}\right].$
(1)

This theory has a continuous $U(1)$ symmetry which shifts $\phi$ by a phase,
so $\phi$ is a natural order parameter for this theory. Notice that the
minimum of the potential is not invariant under this symmetry, so classically
we would expect that the vacuum expectation of $\phi$ to be
$|\left<\phi\right>|=v$. However this may not be the case in the quantum field
theory. To compute the true expectation value of $\phi$ we need to introduce a
symmetry breaking term to the action

$\delta S_{h}=\int d^{d}x~{}h\phi,$ (2)

and then compute the expectation value of $\phi$ in the presence of this term.
The true expectation value of $\phi$ will be given by the limit

$\left<\phi\right>=\lim_{h\rightarrow
0}\left<\phi\right>_{h}=\lim_{h\rightarrow 0}\int D\phi
D\phi^{\dagger}~{}\phi~{}e^{-S-\delta S_{h}}.$ (3)

Though this approach was used to prove the Mermin-Wagner-Coleman theorem for
subsystem symmetries in a lattice theory Batista:2004sc , it not clear how to
define an order parameter for a subsystem symmetry in a continuum field
theory, or how to introduce a symmetry breaking term that only breaks the
desired subsystem symmetry.

A different approach to diagnosing spontaneous symmetry breaking is to compute
the long range correlation functions of charged operators, and see if they
vanish. This approach is used in the context of higher form symmetries to
diagnose if a $p$-form symmetry is spontaneously broken Gaiotto_2015 ;
lake2018higher . For the above example of the complex scalar field $\phi$, one
can compute the two point function

$\left<\phi^{\dagger}(0)\phi(r)\right>.$ (4)

If the field theory is local then the expectation values of a product of
distant operators should factorize, and thus

$\lim_{r\rightarrow\infty}\left<\phi^{\dagger}(0)\phi(r)\right>=\big{<}\phi^{\dagger}\big{>}\big{<}\phi\big{>}=\big{<}\phi\big{>}^{\dagger}\big{<}\phi\big{>}.$
(5)

Hence if $\lim_{r\rightarrow\infty}\left<\phi^{\dagger}(0)\phi(r)\right>=0$
then the symmetry is not spontaneously broken, while if
$\lim_{r\rightarrow\infty}\left<\phi^{\dagger}(0)\phi(r)\right>\neq 0$ then
$\phi$ acquires an expectation value and the symmetry is spontaneously broken.
We note that the local charged operator $\phi$ may need to be renormalized to
give a finite answer. In such a case the UV renormalization needs to be
distinguished from the IR behavior of the 2-point function.

This same method can be used to identify spontaneous symmetry breaking of
subsystem symmetries, though we must ensure that the long distance correlation
functions used to diagnose symmetry breaking is uncharged under the symmetry.
Thus, for a subsystem symmetry, the correlation functions used to identify
symmetry breaking will not be 2-point functions, but rather higher point
correlation functions. We present an example of such a correlation function in
the following subsection.

###  2.1 Symmetry breaking in the XY-plaquette model

As simple example of identifying spontaneous symmetry breaking of subsystem
symmetries, we shall look at continuum field description of the XY plaquette
model Seiberg:2020bhn ; Seiberg:2020wsg ; paramekanti2002ring . In
Seiberg:2020bhn ; Seiberg:2020wsg  it was claimed that no spontaneous symmetry
breaking happens, even in $3+1$ dimensions where spontaneous symmetry breaking
is permitted. We would like to verify this result using our diagnostic for
spontaneous symmetry breaking.

The original XY-plaquette model in $d+1$ dimensions has a continuum field
theory formulation consisting of a compact scalar field $\phi$, and the action
(in Euclidean signature)

$L=\frac{\mu_{0}}{2}\partial_{0}\phi\partial_{0}\phi+\frac{1}{2\mu}\sum_{1\leq
i<j\leq d}\partial_{i}\partial_{j}\phi\partial_{i}\partial_{j}\phi.$ (6)

This model has the global subsystem shift symmetry which acts on the primary
$1+(d-1)$ dimensional planes and is given by

$\phi\rightarrow\phi+\sum_{i=1}^{d}f_{i}(x_{i}),$ (7)

where $f_{i}(x_{i})$ is an arbitrary function of the $i$’th coordinate.

To diagnose if this symmetry is spontaneously broken we should compute a long-
range correlator of charged operators, which for this theory is the operator
$e^{i\phi}$ (as $\phi$ is not a good operator in the continuum field theory
Seiberg:2020bhn ). Furthermore, we need the correlator we compute to be
uncharged under the global subsystem shift symmetry. The simplest correlator
that meets this criteria is the equal time correlator of four $e^{\pm i\phi}$
operators inserted on one of the spatial planes

$\left<e^{i[\phi(0,0,0)-\phi(0,x_{1},0)-\phi(0,0,x_{2})+\phi(0,x_{1},x_{2})]}\right>,$
(8)

where without loss of generality we took the operators to lie in
$x_{1}$–$x_{2}$ plane, and where the remaining spatial coordinates are equal.

As this is a Gaussian theory, we can take the functional integral directly,
and we find that

$\left<e^{i[\phi(0,0,0)-\phi(0,x_{1},0)-\phi(0,0,x_{2})+\phi(0,x_{1},x_{2})]}\right>=e^{-\frac{2}{\mu_{0}}K_{d}(x_{1},x_{2})},$
(9)

where

$\displaystyle K_{d}(x,y)$
$\displaystyle=\int\frac{d\omega}{2\pi}\int\frac{d^{d}k}{(2\pi)^{d}}\frac{(1-\cos(k_{1}x))(1-\cos(k_{2}y))}{\omega^{2}+\frac{1}{\mu\mu_{0}}\sum_{i\neq
j}k_{i}^{2}k_{j}^{2}}$ (10)
$\displaystyle=\frac{\sqrt{\mu\mu_{0}}}{2}\int\frac{d^{d}k}{(2\pi)^{d}}\frac{(1-\cos(k_{1}x))(1-\cos(k_{2}y))}{\sqrt{\sum_{i\neq
j}k_{i}^{2}k_{j}^{2}}}$

Starting with the case $d=2$, we have that

$K_{2}(x,y)=\frac{\sqrt{\mu\mu_{0}}}{8\pi^{2}}\int
d^{2}k\frac{(1-\cos(k_{1}x))(1-\cos(k_{2}y))}{|k_{1}k_{2}|}.$ (11)

This integral is UV divergent, and so must be regulated. As the integral
factorizes, it is natural to introduce independent cutoffs in the $x$ and $y$
directions so that this integral becomes 222We note that the form of (12)
should not depend on the particular regularization scheme. For example on can
compute $K_{2}(x,y)$ using the standard spherical momentum cutoff of
$k_{1}^{2}+k_{2}^{2}\leq\Lambda^{2}$ and find that
$K_{2}(x,y)=\frac{\sqrt{\mu\mu_{0}}}{2\pi^{2}}\left(\gamma+\log(x\Lambda)\right)\left(\gamma+\log(y\Lambda)\right)-\frac{\sqrt{\mu\mu_{0}}}{48}+O(1/\Lambda)$,
which agrees with (12) up to an additive constant.

$\displaystyle K_{2}(x,y)$
$\displaystyle=\frac{\sqrt{\mu\mu_{0}}}{2\pi^{2}}\int_{0}^{\Lambda_{x}}dk_{1}\frac{(1-\cos(k_{1}x))}{k_{1}}\int_{0}^{\Lambda_{y}}dk_{2}\frac{(1-\cos(k_{2}y))}{k_{2}}$
(12)
$\displaystyle=\frac{\sqrt{\mu\mu_{0}}}{2\pi^{2}}\left(\gamma+\log(x\Lambda_{x})+O(1/\Lambda_{x})\right)\left(\gamma+\log(y\Lambda_{y})+O(1/\Lambda_{y})\right).$

This UV divergence cannot be regulated by a renormalization of the operators
(or equivalently by local counter term) due to the presence of the terms
proportional to $\log(x)\log(\Lambda_{y})$ and $\log(y)\log(\Lambda_{x})$. Not
only do these terms make the result unregularizable, their presence indicates
a UV-IR mixing in the charged states of this model, as observed in
Seiberg:2020bhn ; Seiberg:2020wsg . This same correlation function was
computed in the lattice theory in paramekanti2002ring  and
gorantla2021lowenergy , and both found the same result.333One can directly
compare (12) to equation (66) in paramekanti2002ring  and equation (5.33) in
gorantla2021lowenergy .

We can do a similar calculation in $d=3$, where we can take the $k_{3}$
integral by introducing a cutoff, so

$\displaystyle K_{3}(x,y)$
$\displaystyle=\frac{\sqrt{\mu\mu_{0}}}{16\pi^{3}}\int
d^{3}k\frac{(1-\cos(k_{1}x))(1-\cos(k_{2}y))}{\sqrt{k_{1}^{2}k_{2}^{2}+k_{3}^{2}(k_{1}^{2}+k_{2}^{2})}}$
(13) $\displaystyle=\frac{\sqrt{\mu\mu_{0}}}{16\pi^{3}}\int
d^{2}k\frac{(1-\cos(k_{1}x))(1-\cos(k_{2}y))}{\sqrt{k_{1}^{1}+k_{2}^{2}}}\left(2\log(\Lambda_{z})+\log\left(\frac{4(k_{1}^{2}+k_{2}^{2})}{k_{1}^{2}k_{2}^{2}}\right)+O(\Lambda_{z}^{-1})\right).$

As before, this integral will include mixed terms of the form
$f(x,y)\log(\Lambda_{z})$ which cannot be regularized by a local counter term.

In both these instances this diagnostic of spontaneous symmetry breaking
agrees with the results in Seiberg:2020bhn ; Seiberg:2020wsg ;
gorantla2021lowenergy , namely that the symmetry is not spontaneously broken.
Furthermore, these results confirm that the charged states have energies of
order the lattice spacing and so correlation functions involving charged
operators cannot be regularized to give a consistent IR result.

We would like to construct theories with a spontaneously broken continuous
subsystem symmetry, as the IR physics of such theories would be described by
the Goldstone modes of the broken symmetry. Finding such theories would be a
first step to studying the phase transition between the broken and unbroken
phases, a new type of critical phenomena that has yet to be considered in the
literature. In the next two section we will analyze a family of statistical
models inspired by the XY-plaquette model which do indeed have spontaneously
broken subsystem symmetries.

##  3 The classical XY-plaquette model

To realize a theory with spontaneously broken subsystem symmetries, we can try
to start from the lattice model related to the XY Plaquette model in
$d$-dimensions paramekanti2002ring . Consider a square lattice in $d$
dimensions with an angular variable $\phi_{i}$ on each lattice site. In the
original XY model we consider the classical statistical mechanics model where
the energy of each configuration of spins is given by

$E_{\text{XY}}(\phi)=-J\sum_{\left<i,j\right>}\cos(\phi_{i}-\phi_{j}),$ (14)

where the sum runs over all sites $i,j$ which are nearest neighbors on the
lattice. The partition function of this model is given by the path integral

$Z(\beta)=\prod_{i}\int_{0}^{2\pi}d\phi_{i}~{}~{}e^{-\beta
E_{\text{XY}}(\phi)}.$ (15)

As the partition function is given by a Euclidean path integral, we can
consider this model as a quantum mechanical model with a Euclidean action
$S_{E}=E$, while the $d-1$ dimensional Hamiltonian will be related to this
action by a Legendre transform.

We can follow the exact same procedure to construct an XY-plaquette model by
taking the interactions to be on plaquettes rather than between nearest
neighbor sites. This will result in the energy

$E_{\text{XY-plaquette}}(\phi)=-J\sum_{S}\cos(\Delta_{S}\phi),$ (16)

where the sum runs over all square plaquettes $S$ containing the sites
$\\{s_{1},s_{2},s_{3},s_{4}\\}$ in cyclic order,
$\Delta_{S}\phi=\phi_{s_{1}}-\phi_{s_{2}}+\phi_{s_{3}}-\phi_{s_{4}}$, and $J$
is the interaction strength. The partition function of the statistical theory
is described by the Euclidean path integral

$Z(\beta)=\prod_{i}\int_{0}^{2\pi}d\phi_{i}~{}~{}e^{-\beta E_{\text{XY-
plaquette}}(\phi)}.$ (17)

This model is similar to the XY-plaquette model considered in
paramekanti2002ring  and Seiberg:2020bhn , as the plaquette interaction term
is the same, though in our case $E_{\text{XY-plaquette}}(\phi)$ is the
Euclidean action of the quantum model rather than the Hamiltonian, and it does
not contain a kinetic term.

This XY-plaquette model has the subsystem symmetry of shifting all
$\phi_{i}$’s on any $d-1$ dimensional hyper-plane in the lattice by a constant
angle. Similar to the original XY model, we would expect this model to be in a
broken phase at low energy if the dimension allows for such a phase to exist.
The continuum limit of this broken phase should correspond to taking the
continuum limit of the of the Euclidean action on the lattice. This was done
carefully in Seiberg:2020bhn , and the resulting model is a compact scalar
field $\phi$ in $d$ dimensions with the Lagrangian

$L=\frac{1}{2g}\sum_{1\leq i<j\leq
d}\partial^{i}\partial^{j}\phi\partial_{i}\partial_{j}\phi.$ (18)

This Lagrangian has the global dipole symmetry of shifting $\phi$ by any
arbitrary function of a single variable

$\phi\rightarrow\phi+\sum_{i=1}^{d}f_{i}(x_{i}).$ (19)

As $\phi$ is $2\pi$ periodic, these functions are defined up to an integer
multiple of $2\pi$. Thus neither $\phi$ nor $\partial_{i}\phi$ are good
operators in the continuum limit, while $\partial_{i}\partial_{j}\phi$ is a
good operator so long as $i\neq j$, as is $e^{i\phi}$, similar to the analysis
of Seiberg:2020bhn . Furthermore, we will only focus on field configurations
that have a finite action in the continuum limit. This excludes any winding
modes that have infinite action (due to a computation that is identical to the
computation done in section 4.2 of Seiberg:2020bhn ).444In this model the
momentum modes have finite action, as will be shown bellow, and so we can
consistently exclude the winding modes from the continuum field theory. In
Seiberg:2020bhn  both momentum and winding modes had infinite action in the
continuum limit, and so one needs to keep both.

The spontaneous breaking of this dipole symmetry is what we are interested in
analyzing. As this subsystem symmetry acts on $d-1$ dimensional hyper-planes,
the minimal dimension we can hope to have a spontaneously broken symmetry is
$d=4$, and indeed we will see that this is the case. However before we analyze
the breaking of the dipole symmetry we present an analysis of the spectrum of
the model. We note that a further study the correlation functions of uncharged
operators is presented in appendix A.

###  3.1 Spectrum

To understand the spectrum of the model, we start from the equation of motion
for $\phi$

$\sum_{1\leq i<j\leq d}\partial_{i}^{2}\partial_{j}^{2}\phi=0.$ (20)

Then we can expand $\phi$ in terms of Fourier modes, which leads to the
dispersion relation

$\sum_{1\leq i<j\leq d}k_{i}^{2}k_{j}^{2}=0.$ (21)

If we consider the theory in Lorentzian signature, that is take one of the $d$
spatial coordinates and Wick rotate it, then the dispersion relation becomes

$\omega^{2}=\frac{\sum_{i\neq j}k_{i}^{2}k_{i}^{2}}{k^{2}}.$ (22)

We note that this dispersion relation is continuous around $k=0$, and indeed
$\omega\rightarrow 0$ as $|k|\rightarrow 0$ from any direction.

Next we quantize the Lorenzian theory on a $d-1$-torus of equal length $\ell$,
so that the momenta are quantized $k_{i}=2\pi\frac{n_{i}}{\ell}$. In terms of
the momentum modes

$\phi_{\mathbf{n}}=\frac{1}{(2\pi)^{d-1}}\int d^{d-1}x~{}\phi(x)e^{2\pi
i\frac{\mathbf{n}\cdot\mathbf{x}}{\ell}}$ (23)

the Lagrangian becomes

$L=\frac{\ell^{d-3}(2\pi)^{2}}{2g}\sum_{\mathbf{n}\in\mathbb{Z}^{3}}\left[|\mathbf{n}|^{2}\partial_{0}\phi_{-\mathbf{n}}\partial_{0}\phi_{\mathbf{n}}-\frac{(2\pi)^{2}}{\ell^{2}}\sum_{i<j}n_{i}^{2}n_{j}^{2}\phi_{-\mathbf{n}}\phi_{\mathbf{n}}\right].$
(24)

Each of these modes for generic values on $\mathbf{n}$ behaves as a simple
harmonic oscillator with ground state energy

$E_{\mathbf{n}}=\frac{1}{2}\omega_{\mathbf{n}}=\frac{\pi}{\ell}\left(\frac{\sum_{i\neq
j}n_{i}^{2}n_{i}^{2}}{\mathbf{n}^{2}}\right)^{1/2},$ (25)

where $\omega$ given by the dispersion relation (22). The rest of the states
are built via the generic Fock space construction on these ground states. If,
however, all but one of the integers in $\mathbf{n}$ are zero then
$\omega_{\mathbf{n}}=0$ and these fields have no harmonic potential. We may
need to quantize these modes more carefully, as they can have large momenta in
the remaining direction.

To quantize the Lagrangian for the modes where all but one of the
$\mathbf{n}$’s are zero, we follow Seiberg:2020wsg  and write these modes in
position space by writing the field as

$\phi(x_{0},x_{i},\ldots,x_{d-1})=\phi_{1}(x_{0},x_{1})+\phi_{2}(x_{0},x_{2})+\ldots+\phi_{d-1}(x_{0},x_{d-1})+\ldots$
(26)

Then the Lagrangian for these modes is

$L=\frac{1}{2g}\sum_{i=1}^{d-1}\left(\partial_{0}\partial_{i}\phi_{i}\right)^{2}.$
(27)

We see that unlike in Seiberg:2020wsg  these modes are decoupled, and the
“gauge” symmetry that connects these modes is just a subset of the global
symmetry of the Lagrangian $\phi\rightarrow\phi+f(t)$. Overall these modes
behave just like the zero modes of a standard periodic free scalar field, and
do not acquire a mass of the order of the lattice spacing, which is very
different than the theories considered in Seiberg:2020bhn ; Seiberg:2020wsg .

This analysis shows that the even for the modes where all but one of the
$\mathbf{n}$’s are zero the Lagrangian in (24) is the correct Lagrangian, and
the Hamiltonian for these modes in the $i$’th direction is

$H_{i}=\frac{\ell^{d-3}(2\pi)^{2}}{2g}\sum_{n\neq
0}\frac{1}{n^{2}}\pi_{i,-n},\pi_{i,n}.$ (28)

The conjugate momenta of these modes, denoted by $\pi_{i,n}$, are the
generators of the global symmetry $\phi\rightarrow\phi+f_{i}(x_{i})$. These
conjugate momenta do indeed commute with the Hamiltonian, as expected for a
global symmetry, though their spectrum is continuous. Choosing a ground state
for the system is the same as choosing the eigenvalues of these modes, and so
the $\pi_{i,n}$’s shift us from one ground state to another, in analogy to
typical Goldstone modes. Thus this picture looks like a spontaneously broken
symmetry, with the $\phi$ being an analogous Goldstone boson.

There is also the special mode $\mathbf{n}=0$, for which the Lagrangian
vanishes and so is not a physical mode. This is very different than in
Seiberg:2020wsg  where this mode is the unique ground state of the system. The
reason this mode disappears from the theory is the time translation symmetry
$\phi\rightarrow\phi+f(t)$, which is generated by this zero mode.

Finally we note is that the resulting Hamiltonian is non-local in space for
all of the modes. Indeed for the continuum theory the conjugate variable to
$\phi$ is

$\pi=\frac{\partial
L}{\partial(\partial_{0}\phi)}=\frac{1}{g}\partial_{0}\nabla^{2}\phi,$ (29)

so the Hamiltonian of the system is

$H=\int d^{d-1}x~{}\left[\frac{1}{2g}\pi\nabla^{-2}\pi+\frac{1}{2g}\sum_{1\leq
i<j\leq d-1}\left(\partial_{i}\partial_{j}\phi\right)^{2}\right],$ (30)

which has a non-local kinetic term. This Hamiltonian still has a fine
expansion in terms of Fourier modes, as seen above, and so can be quantized in
the standard way.

It is interesting that a local Lagrangian gave rise to a non-local
Hamiltonian, which raises the question which is more fundamental. From an
effective field theory prospective, one writes down the most general local
Lagrangian that obeys the desired symmetries. Then one keeps only the
renormalizable terms, which are found by computing the scaling dimensions.
Based on this philosophy a local Lagrangian is more fundamental to
understanding the IR physics, as it contains the relevant information about
the symmetries and scaling dimensions. The model we presented fits nicely into
this picture, as it is the simplest Lagrangian that obeys the desired
subsystem symmetry, and indeed it seems that the resulting quantum field
theory is local and well defined.

There is however the alternate bottom up approach where one starts from a
local UV Hamiltonian, say on a lattice, and then tries to infer the physics at
large distances. In this case one assumes that the Hamiltonian is local, and
so it does not seem that the classical XY model can describe the IR physics
that emerge from such a construction. This raises many interesting questions
relating to the existence of phases of matter that fit into an effective field
theory picture, but that cannot be described by a local lattice Hamiltonian.

There is also the possibility that certain phases of matter may arise from a
local lattice Hamiltonian that has a non-local Lagrangian. For example one can
introduce a nearest neighbor interaction of the conjugate momenta variables to
a lattice Hamiltonian, which would keep the Hamiltonian local but result in a
non-local Lagrangian. Such systems have been considered in the condensed
matter literature, for example Wu_2021  studied similar Hamiltonians with
subsystem symmetries, though it is not clear to us how such models fit into
the canonical picture of Wilsonian quantum field theory.

###  3.2 Spontaneous subsystem symmetry breaking

To identify if the subsystem symmetry (19) is spontaneously broken, we will
look at long range correlation functions of the charged local operator
$e^{i\phi}$.555We explore the correlation functions of uncharged operators and
their consequences in appendix A. The 1, 2, and 3-point functions of these
operators must vanish by the global dipole symmetry. Hence the simplest
nonzero correlation function of such operators is the 4-point rectangle
function, similar to what was considered in (8). This correlator is

$\left<e^{i[\phi(0,0)-\phi(x_{1},0)-\phi(0,x_{2})+\phi(x_{1},x_{2})]}\right>,$
(31)

where we have chosen the rectangle to lie in the $(x_{1},x_{2})$ plane, and
suppressed the dependence on the remaining $d-2$ coordinates.

Taking the Gaussian functional integral we can write this expectation value as

$\left<e^{i[\phi(0,0)-\phi(x_{1},0)-\phi(0,x_{2})+\phi(x_{1},x_{2})]}\right>=e^{-2gI_{d}(x_{1},x_{2})},$
(32)

where

$I_{d}(x,y)=\int\frac{d^{d}k}{(2\pi)^{d}}\frac{(1-\cos(k_{1}x))(1-\cos(k_{2}y))}{\sum_{i\neq
j}k_{i}^{2}k_{j}^{2}}.$ (33)

Notice that this integral is well defined in the IR, as around $k_{1},k_{2}\ll
1$ this integral behave as $(k_{1}^{2}k_{2}^{2})/(\sum_{i\neq
j}k_{i}^{2}k_{j}^{2})$, which is finite around zero for a fixed angle.

For $d=2$ we can take this integral analytically, and we find

$I_{2}(x,y)=\frac{|xy|}{4}.$ (34)

This grows as we take $x$ and $y$ to be large, hence the full correlator (31)
vanishes, and thus the subsystem symmetry is unbroken.

For $d=3$ this integral evaluates to

$\displaystyle I_{3}(x,y)$
$\displaystyle=\frac{1}{4\pi}\left[|y|\log\left(\frac{|x|+\sqrt{x^{2}+y^{2}}}{|y|}\right)+|x|\log\left(\frac{|y|+\sqrt{x^{2}+y^{2}}}{|x|}\right)\right]$
(35)
$\displaystyle=\frac{1}{4\pi}|xy|^{1/2}\left\\{\alpha^{-1/2}\log\left(\alpha+\sqrt{1+\alpha^{2}}\right)+\alpha^{1/2}\log\left(\frac{1+\sqrt{1+\alpha^{2}}}{\alpha}\right)\right\\},$

where $\alpha\equiv\frac{|x|}{|y|}$. Again the integral grows as $x$ or $y$
(or both) are taken to be large, implying that the subsystem symmetry is still
unbroken. This is expected as the minimal dimension allowed for spontaneous
symmetry breaking in this model is $d=4$.

For $d=4$ we can write the rectangular function as

$\displaystyle I_{4}(x,y)$ $\displaystyle=\frac{1}{8\pi^{3}}\int
dk_{1}~{}dk_{2}~{}\frac{(1-\cos(k_{1}x))(1-\cos(k_{2}y))}{k_{1}^{2}+k_{2}^{2}}K\left(1-\frac{k_{1}^{2}k_{2}^{2}}{(k_{1}^{2}+k_{2}^{2})^{2}}\right),$
(36)

where $K(m)$ is the complete elliptic integral of the first kind. Notice that
$I_{4}(x,y)$ is only a function of the ratio $x/y$, and so we expect it to be
a constant for large rectangles with the ratio $x/y$ fixed. However this
integral is logarithmically divergent in the UV, and this UV divergence neads
to be regularized. Moving to polar coordinates and introducing a UV cuttoff,
we can take the radial integral to find

$\displaystyle I_{4}(x,y)$
$\displaystyle=\frac{1}{8\pi^{3}}\int_{0}^{\Lambda}dr\int_{0}^{2\pi}d\theta~{}\frac{(1-\cos(xr\cos\theta))(1-\cos(yr\sin\theta))}{r}K\left(1-\cos^{2}\theta\sin^{2}\theta\right)$
(37)
$\displaystyle=\frac{1}{8\pi^{3}}\int_{0}^{2\pi}d\theta~{}K\left(1-\cos^{2}\theta\sin^{2}\theta\right)\left[\log\Lambda+\gamma+\frac{1}{2}\log\left(\frac{x^{2}y^{2}\cos^{2}\theta\sin^{2}\theta}{\left|x^{2}\cos^{2}\theta-y^{2}\sin^{2}\theta\right|}\right)\right]+O(\Lambda^{-1}).$

Notice that the angular integral multiplying the $\log\Lambda$ term is finite,
hence we can renormalize the operator by a local counterterm to cancel the
logarithmic UV divergence, leaving the remaining integral finite and cutoff
independent. Canceling the UV cutoff does however make the integral
logarithmically diverge as the area of the rectangles is taken to be large
while fixing the ratio of the two sides. However this term approaches a
constant if only one side of the rectangle is taken to be large while the
other is kept fixed (notice that this is not the case for $I_{2}$ and $I_{3}$
which both diverge even in this limit, though slower than the divergence for
fixed ratio). This implies that the subsystem symmetry is indeed spontaneously
broken for $d=4$.

We expect a similar behavior in dimensions $d>4$, where an operator
renormalization can take care of the UV divergence, while no IR divergent will
exist at all. This is in alignment with the previous discussions of
spontaneously breaking subsystem symmetries based on Batista:2004sc . For
completeness we present plots of $I_{3}(x,y)$ and $I_{4}(x,y)$ in figures 1
and 2.

Figure 1: A surface plot of $I_{3}(x,y)$, as well as a plot of $I_{3}(x,y)$
for fixed values of $y$. We also compare $I_{3}(x,y)$ at fixed $y$ to the
asymptotic behavior at large $x$ denoted by the dotted lines. The asymptotic
behavior is
$I_{3}(x,y)\approx\frac{y}{4\pi}\left(\log(2x/y)+1+O(1/x^{2})\right)$.  Figure
2: A surface plot of $I_{4}(x,y)$, as well as a plot of $I_{4}(x,y)$ for fixed
values of $y$. We also compare $I_{4}(x,y)$ at fixed $y$ to the asymptotic
value at infinite $x$ denoted by the dotted lines.

##  4 Other classical XY models with subsystem symmetries

The classical XY-plaquette model can be generalized to a family of models in
$d$ dimensions that has a subsystem symmetry which acts on $d-m+1$ dimensional
hyperplanes. As before we start with a square lattice in $d$ dimensions with
an angular variable $\phi_{i}$ on each lattice site, though now the energy of
each configuration is given by

$E(\phi)=-J\sum_{C_{m}}\cos(\Delta_{C_{m}}\phi).$ (38)

This the sum runs over all hyper-cubes $C_{m}$ containing the sites
$\\{c_{s+\sum_{j=1}^{m}i_{j}\mathbf{e}_{j}}\\}$, where $s$ is any lattice
site, $i_{j}\in\\{0,1\\}$, and $\mathbf{e}_{m}$ are $m$ orthogonal unit
vectors on the lattice. The classical XY model is the case $m=1$, while the
classical XY-plaquette model is the case $m=2$.

The next class of models is the case $m=3$, which is a classical version of
the XY-cube model considered in You_2020 ; Gorantla:2020xap . We present a
brief analysis of this model below, before considering the general $m$ case.

###  4.1 The classical XY-cube model

The classical XY-cube model has the subsystem symmetry of shifting all
$\phi_{i}$’s on any $d-2$ dimensional hyper-plane in the lattice by a constant
angle. Similar to the previous model, we would expect this model to be in a
broken phase at low energy if the dimension allows for such a phase to exist,
which for this model is $d>4$. The continuum limit of this broken phase should
correspond to taking the continuum limit of the of the Euclidean action on the
lattice. This was done carefully in Gorantla:2020xap , and the resulting model
is a compact scalar field $\phi$ in $d$ dimensions with the Lagrangian

$L=\frac{1}{2g}\sum_{1\leq i<j<k\leq
d}\partial^{i}\partial^{j}\partial^{k}\phi\partial_{i}\partial_{j}\partial_{k}\phi.$
(39)

This Lagrangian has the global dipole symmetry of shifting $\phi$ by any
arbitrary function of any two variables

$\phi\rightarrow\phi+\sum_{1\leq i<j\leq d}f_{i,j}(x_{i},x_{j}).$ (40)

As $\phi$ is $2\pi$ periodic, these functions are defined up to an integer
multiple of $2\pi$. Thus neither $\phi$ nor $\partial_{i}\phi$ or
$\partial_{i}\partial_{j}\phi$ are good operators in the continuum limit,
while $\partial_{i}\partial_{j}\partial_{k}\phi$ is a good operator so long as
$i\neq j\neq k$, as is $e^{i\phi}$, similar to the analysis of
Gorantla:2020xap . As before, we will only focus on field configurations that
have a finite action in the continuum limit, excluding any winding modes that
have infinite action.

To find the spectrum of the theory, we can expand $\phi$ in terms of Fourier
modes, which leads to the dispersion relation

$\sum_{1\leq i<j<l\leq d}k_{i}^{2}k_{j}^{2}k_{l}^{2}=0.$ (41)

If we consider the theory in Lorentzian signature, that is take one of the $d$
spatial coordinates and Wick rotate it, then the dispersion relation becomes

$\omega^{2}=\frac{\sum_{i<j<l}k_{i}^{2}k_{k}^{2}k_{l}^{2}}{\sum_{i<j}k_{i}^{2}k_{j}^{2}}.$
(42)

As before, this dispersion relation is continuous around $k=0$.

We can quantize the Lorenzian theory on a $d-1$-torus of equal length $\ell$,
so that the momenta are quantized $k_{i}=2\pi\frac{n_{i}}{\ell}$. In terms of
the momentum modes given in (23), the Lagrangian becomes

$L=\frac{\ell^{d-5}(2\pi)^{4}}{2g}\sum_{\mathbf{n}\in\mathbb{Z}^{3}}\left[\left(\sum_{i<j}n_{i}^{2}n_{j}^{2}\right)\partial_{0}\phi_{-\mathbf{n}}\partial_{0}\phi_{\mathbf{n}}-\frac{(2\pi)^{2}}{\ell^{2}}\sum_{i<j<k}n_{i}^{2}n_{j}^{2}n_{k}^{2}\phi_{-\mathbf{n}}\phi_{\mathbf{n}}\right].$
(43)

Each of these modes for generic values on $\mathbf{n}$ behaves as a simple
harmonic oscillator with ground state energy

$E_{\mathbf{n}}=\frac{1}{2}\omega_{\mathbf{n}}=\frac{\pi}{\ell}\left(\frac{\sum_{i<j<k}n_{i}^{2}n_{j}^{2}n_{k}^{2}}{\sum_{i<j}n_{i}^{2}n_{j}^{2}}\right)^{1/2},$
(44)

with $\omega$ given by the dispersion relation (42). The rest of the states
are built via the generic Fock space construction on these ground states. If,
however, all but two of the integers in $\mathbf{n}$ are zero then $\omega=0$
and these modes have no harmonic potential. As before, these modes behave just
like the zero modes of a standard periodic free scalar field, and do not
acquire a mass of the order of the lattice spacing, which is very different
than the theories considered in Seiberg:2020bhn ; Seiberg:2020wsg ;
Gorantla:2020xap .

This analysis shows that, even for the modes where all but two of the
$\mathbf{n}$’s are zero, the Lagrangian in (43) is correct. The momenta
conjugate to these modes, denoted by $\pi_{i,j,n,m}$, are the generators of
the global symmetries $\phi\rightarrow\phi+f_{ij}(x_{i},x_{j})$. These
conjugate momenta do indeed commute with the Hamiltonian, as expected for a
global symmetry. Choosing a ground state for the system is the same as
choosing the eigenvalues of these operators, and so the $\pi_{i,j,n,m}$’s
shift us from one ground state to another, just as the zero modes for the
usual Goldstone boson. The modes corresponding to where all but one of the
$n$’s are zero now have vanishing action, and correspond to the symmetries
involving the time direction $\phi\rightarrow\phi+f_{0j}(x_{0},x_{j})$.
Finally we note that the resulting Hamiltonian is non-local in space, and the
discussion regarding this non-locality from section 3.1 applies to this model
as well.

As in the XY-Plaquette model, to check for which dimensions the symmetry is
spontaneously broken we need to compute correlation functions of charged
operators. In this model the simplest invariant correlation function is the
cubic correlator of eight $e^{i\phi}$ given by

$\left<\exp\left[i\sum_{i,j,k\in\\{0,1\\}}(-1)^{i+j+k}\phi(ix,jy,kz)\right]\right>.$
(45)

As this theory is Gaussian, we can take the path integral directly resulting
in

$\left<\exp\left[i\sum_{i,j,k\in\\{0,1\\}}(-1)^{i+j+k}\phi(ix,jy,kz)\right]\right>=e^{-4gJ_{d}(x,y,z)}$
(46)

where

$J_{d}(x,y,z)=\int\frac{d^{d}k}{(2\pi)^{d}}\frac{(1-\cos(k_{1}x))(1-\cos(k_{2}y))(1-\cos(k_{3}z))}{\sum_{i<j<l}k_{i}^{2}k_{j}^{2}k_{l}^{2}}.$
(47)

For $d=3$, the lowest dimension for which the model makes sense, we have that

$J_{3}(x,y,z)=\frac{1}{8}|xyz|$ (48)

As before, this implies that the symmetry is not broken as the long range
correlation functions grows.

For $d=4$ this integral is

$J_{4}(x,y,z)=\int\frac{d^{3}k}{16\pi^{3}}\frac{(1-\cos(k_{1}x))(1-\cos(k_{2}y))(1-\cos(k_{3}z))}{|k_{1}k_{2}k_{3}|\sqrt{k_{1}^{2}k_{2}^{2}+k_{2}^{2}k_{3}^{2}+k_{1}^{2}k_{3}^{2}}}.$
(49)

This integral is completely convergent for any finite values of $x,y,z$,
though a closed form expression for the integrand is beyond our reach. We can
however take various large distance limits of this function. Just via scaling
arguments it is clear that if we take $x,y$ and $z$ to be large while keeping
their ratios fixed then

$J_{4}(x\rightarrow\infty,y\rightarrow\infty,z\rightarrow\infty)=|xyz|^{2/3}f(x/y,x/z),$
(50)

where

$f(a,b)=\int\frac{d^{3}k}{16\pi^{3}}\frac{(1-\cos(k_{1}))(1-\cos(k_{2}))(1-\cos(k_{3}))}{|k_{1}k_{2}k_{3}|\sqrt{a^{2/3}b^{-4/3}k_{1}^{2}k_{2}^{2}+a^{2/3}b^{2/3}k_{2}^{2}k_{3}^{2}+a^{-4/3}b^{2/3}k_{1}^{2}k_{3}^{2}}}.$
(51)

If we take only $x$ and $y$ to be large while keeping $z$ and the ratio $x/y$
fixed then

$\displaystyle J_{4}(x\rightarrow\infty,y\rightarrow\infty,z)$
$\displaystyle=\frac{1}{4}|xy|^{1/2}|z|\int\frac{d^{2}k}{(2\pi)^{2}}\frac{(1-\cos
k_{1})(1-\cos k_{2})}{|k_{1}k_{2}|\sqrt{(y/x)k_{1}^{2}+(x/y)k_{2}^{2}}}$ (52)
$\displaystyle=\frac{|z||xy|^{1/2}}{8\pi}\left\\{\alpha^{-1/2}\log\left(\alpha+\sqrt{1+\alpha^{2}}\right)+\alpha^{1/2}\log\left(\frac{1+\sqrt{1+\alpha^{2}}}{\alpha}\right)\right\\},$

where $\alpha\equiv|x/y|$.

Finally, taking $x\rightarrow\infty$ while keeping $y$ and $z$ fixed we see
that the integral scales as

$J_{4}(x\rightarrow\infty,y,z)=\frac{1}{8\pi}|yz|\big{(}\log|x|+O(1)\big{)},$
(53)

so even in this limit the cube correlation function diverges. As this symmetry
acts on $d-2$ dimensional hyper-planes, it cannot be spontaneously broken in
dimensions $d\leq 4$, so this analysis does in fact make sense.

For this symmetry to be spontaneously broken we need to be in at least $d=5$,
where the integral evaluates to

$\displaystyle J_{5}(x,y,z)$
$\displaystyle=\int\frac{d^{4}k}{2(2\pi)^{4}}\frac{(1-\cos(k_{1}x))(1-\cos(k_{2}y))(1-\cos(k_{3}z))}{\sqrt{\left(\sum_{i<j<l}k_{i}^{2}k_{j}^{2}k_{l}^{2}\right)\left(\sum_{i<j}k_{i}^{2}k_{j}^{2}\right)}}$
(54)
$\displaystyle=\int\frac{d^{3}k}{(2\pi)^{4}}\frac{(1-\cos(k_{1}x))(1-\cos(k_{2}y))(1-\cos(k_{3}z))}{|k_{1}k_{2}k_{3}|\sqrt{k_{1}^{2}+k_{2}^{2}+k_{3}^{2}}}K\left(1-\frac{\left(k_{1}^{2}k_{2}^{2}+k_{1}^{2}k_{3}^{2}+k_{2}^{2}k_{3}^{2}\right)^{2}}{k_{1}^{2}k_{2}^{2}k_{3}^{2}\left(k_{1}^{2}+k_{2}^{2}+k_{3}^{2}\right)}\right).$

This integral is still UV finite, and converges for any values of $x,y$ and
$z$, though it has no known closed form expression. We can however do the same
limit analysis to how it behaves when we take the cube to be large. If we take
$x\rightarrow\infty$ while keeping $y$ and $z$ fixed then the integral
approximates as 666We can write $J_{5}=\int dk_{1}(1-\cos(k_{1}))f(k_{1})$,
and as $f(k_{1})$ is an integrable function its Fourier transform decays at
infinity. Thus the integral $\int dk_{1}\cos(k_{1}x)f(k_{1})$ vanishes as
$x\rightarrow\infty$ and we are left with the remaining finite integral.

$\displaystyle J_{5}(x\rightarrow\infty,y,z)=2\pi\int
d^{3}k\frac{(1-\cos(k_{2}y))(1-\cos(k_{3}z))}{|k_{1}k_{2}k_{3}|\sqrt{k_{1}^{2}+k_{2}^{2}+k_{3}^{2}}}K\left(1-\frac{\left(k_{1}^{2}k_{2}^{2}+k_{1}^{2}k_{3}^{2}+k_{2}^{2}k_{3}^{2}\right)^{2}}{k_{1}^{2}k_{2}^{2}k_{3}^{2}\left(k_{1}^{2}+k_{2}^{2}+k_{3}^{2}\right)}\right),$
(55)

which is a constant value independent of $x$. Hence this symmetry is indeed
broken in $d=5$, as expected.

###  4.2 Classical XY-hypercube models

We now turn our attention back to the general $m$ case, the classical XY-
hypercube model of degree $m$, with the energy given by (38). This model has a
subsystem symmetry which acts on $d-m+1$ dimensional hyperplanes. The
continuum field theory of these models consists of a periodic scalar field
$\phi$ with the Lagrangian

$L=\frac{1}{2g}\sum_{i_{1}<i_{2}<\ldots<i_{m}=1}^{d}\partial^{i_{1}}\partial^{i_{2}}\ldots\partial^{i_{m}}\phi~{}\partial^{i_{1}}\partial^{i_{2}}\dots\partial^{i_{m}}\phi.$
(56)

This model has the subsystem shift symmetry given by

$\phi\rightarrow\phi+f_{i_{1},i_{2},\ldots,i_{m-1}}(x_{i_{1}},x_{i_{2}},\ldots,x_{i_{m-1}})$
(57)

which acts on $d-m+1$ dimensional hyperplanes.

The spectrum of the model with arbitrary $m$ has many of the same features as
the spectrum of the plaquette and cubic models. The dispersion relation in
Lorentzian signature will be given by

$\omega^{2}=\frac{\sum_{i_{1}<\ldots<i_{m}}k_{i_{1}}^{2}\ldots
k_{i_{m}}^{2}}{\sum_{i_{1}<\ldots<i_{m-1}}k_{i_{1}}^{2}\ldots
k_{i_{m-1}}^{2}},$ (58)

The quntization of the theory on a torus follows the same procedure as was
done for the plaquette and cubic models. Generic modes will behave like
Harmonic oscillators with ground state energy $E=\frac{1}{2}\omega$, where
$\omega$ is given by the dispersion relation (58). Modes where all but $m-1$
of the $k$’s have zero harmonic potential but a non-zero kinetic term, and so
are the zero modes of the theory. These modes generate the subset of shift
symmetries in (57) that do not involve the time direction. Similarly, the
modes where all but $m-2$ of the $k$’s are zero have a vanishing Lagrangian
and generate the symmetries that do involve the time direction.

The dimensions for which the subsystem symmetry (57) can be spontaneously
broken is $d>1+m$ (by the generalized Merman-Wagner-Coleman theorem
Batista:2004sc ). We have shown that the symmetry is indeed spontaneously
broken when $d=2+m$ for $m\leq 3$777We considered $m=2,3$ in this paper, while
the case $m=1$ is just a free compact scalar field which is known to exhibit
spontaneous symmetry breaking in $d=3$., and expect this trend to continue for
arbitrary $m$.

To diagnose whether the subsystem symmetry (57) is spontaneously broken, one
would need to compute a correlation function of $e^{\pm i\phi}$ operators
placed on the corners of a $m$-dimensional hypercube, as this is the simplest
correlation function of charged operators that is uncharged under the
subsystem symmetry. Calling this correlation function
$F_{m}(x_{1},x_{2},\ldots,x_{m})$, one can preform the Gaussian path integral
to find that

$F_{m}(x_{1},x_{2},\ldots,x_{m})=e^{-2^{m-1}g\;I_{m,d}},$ (59)

where

$I_{m,d}(x_{1},x_{2},\ldots,x_{m})=\int\frac{d^{d}k}{(2\pi)^{d}}\frac{\prod_{i=1}^{m}(1-\cos(k_{i}x_{i}))}{\sum_{i_{1}<\ldots<i_{m}}k_{i_{1}}^{2}\ldots
k_{i_{m}}^{2}}.$ (60)

The behavior of $I_{m,d}(x_{1},x_{2},\ldots,x_{m})$ at large separations
determines the spontaneous breaking of the subsystem symmetry; if $I_{m,d}$
grows at large separation then the symmetry is unbroken while if it approaches
a constant then the symmetry is spontaneously broken. We have shown that
$I_{2,4}$ and $I_{3,5}$ approach a constant when one of the coordinates is
taken to be large, and expect $I_{m,m+2}$ to behave the same way for arbitrary
$m$.

We can actually compute this limit directly from (60). Notice that $I_{m,m}$
factorizes, and thus we can directly compute

$I_{m,m}(x_{1},x_{2},\ldots,x_{m})=2^{-m}|x_{1}x_{2}\ldots x_{m}|.$ (61)

Hence the symmetry is not spontaneously broken for $d=m$.

For $d=m+1$ we can take the limit $x_{1}\rightarrow\infty$ while keeping the
rest of the variables fixed,888The simplest way to do this calculation is to
first take both $x_{1}$ and $x_{2}$ to be large, when the limit reduces to an
integral similar to (52). Then one can take the limit $x_{1}\rightarrow\infty$
while keeping $x_{2}$ fixed and obtain (62). and we find that

$I_{m,m+1}(x_{1}\rightarrow\infty,x_{2},\ldots,x_{m})=2^{-m}\pi^{-1}|x_{2}\ldots
x_{m}|\log(x_{1}),$ (62)

again implying the symmetry is not spontaneously broken.

For $d=m+2$ the integral $I_{m,m+2}$ is UV finite assuming $m>2$, so we can
take the same limit $x_{1}\rightarrow\infty$ while keeping the rest of the
variables fixed, resulting in the finite value 999This limit is obtained in
the same way as (55).

$\displaystyle I_{m,m+2}(x_{1}\rightarrow\infty,x_{2},\ldots,x_{m})=\int$
$\displaystyle\frac{d^{m}k}{(2\pi)^{m+1}}\frac{\prod_{i=2}^{m}(1-\cos(k_{i}x_{i}))}{\prod_{i=1}^{m}|k_{i}|\left(\sum_{i_{1}<\ldots<i_{m-2}}k_{i_{1}}^{2}\ldots
k_{i_{m-2}}^{2}\right)^{1/2}}\times$ (63) $\displaystyle\qquad\times
K\left(1-\frac{\left(\sum_{i_{1}<\ldots<i_{m-1}}k_{i_{1}}^{2}\ldots
k_{i_{m-1}}^{2}\right)^{2}}{\prod_{i=1}^{m}k_{i}^{2}\left(\sum_{i_{1}<\ldots<i_{m-2}}k_{i_{1}}^{2}\ldots
k_{i_{m-2}}^{2}\right)}\right).$

As in this limit $I_{m,m+2}$ approaches a constant (i.e. is independent of
$x_{1}$) the symmetry is indeed spontaneously broken. Hence the subsystem
symmetry (57) is spontaneously broken when $d=m+2$ for arbitrary $m$, but not
when $d\leq m+1$. This is in agreement with the expectations based on the
generalized Merman-Wagner-Coleman theorem Batista:2004sc .

##  5 Discussion and future directions

Our main goal in this paper was to discuss spontaneous breaking of subsystem
symmetries, and to construct well behaved continuum quantum field theories
that model the spontaneously broken phase. Though we successfully built field
theories that exhibit a spontaneous broken subsystem symmetry, there are still
many open questions relating to these models, and to quantum field theories
with subsystem symmetries in general. In this section we discuss in further
detail two general areas of study that we find most interesting. The first is
to study the unbroken phase and the phase transition in the classical XY-
plaquette model we constructed in section 3. The second is to analyze
generalizations and interacting theories built upon the classical XY-plaquette
model.

There are however many other interesting aspects to keep investigating with
regards to the models we constructed. One natural question to ask is whether
the subsystem symmetries are still spontaneously broken at finite temperature,
assuming such breaking is allowed by the number of dimensions. As our models
have no free parameters, the symmetry is either broken or unbroken at all
finite temperatures. To understand which is the case, we need to compute the
same 4-point correlation functions of exponential operators we considered when
taking one of the spatial dimensions to be compact. Though one can evaluate
the sum over the Matsubara frequencies in the compact direction, the resulting
integrals over the momenta in the remaining directions are unwieldy, and must
be UV regularized for the classical XY-plaquette model. Hence we decided to
leave this problem for a future endeavour.

Another interesting observation is that the Lagrangian (18) in $d=2$ is
similar to the Lagrangian of a zero velocity chiral boson,
$\mathcal{L}=\alpha\partial_{0}\phi\partial_{1}\phi$. The chiral boson also
possesses the same subsystem symmetry (19), which is not spontaneously broken
in $d=2$. One can try to generalize the chiral Lagrangian to higher
dimensions, namely consider the Lagrangian $\mathcal{L}=\alpha\sum_{i\neq
j}\partial_{i}\phi\partial_{j}\phi$ which also possesses the same subsystem
symmetry. However this Lagrangian breaks the reflection symmetry of the
lattice, and also has a Hamiltonian that is unbounded from below. One can
overcome some of the difficulties in quantizing the $d=2$ theory using the
conformal symmetry, and perhaps in higher dimensions there is a different
conformal symmetry which emerges, similar to the ones studied in Karch_2021 ,
which may be worth studying.

###  5.1 The unbroken phase and phase transition

Up to now we have focused on the symmetry broken phases of the classical XY
models we considered. However, just like the standard XY model, the lattice
models have a phase where the symmetry is unbroken, and a phase transition
between the two phases. To model the phase transition and the unbroken phase,
we would like to construct a Lagrangian for a complex scalar field $\Phi$ that
has the desired $U(1)$ subsystem symmetry. Then tuning the potential of $\Phi$
in this Lagrangian would (at least classically) give $\Phi$ a vacuum
expectation value, and thus take us between the spontaneously broken and
unbroken phases.

This Lagrangian can be built in a similar manner to the Lagrangian’s with
global dipole symmetries constructed in Pretko_2018 ; Gromov2019 . Focusing on
the classical XY-plaquette model, the lowest order terms in this effective
Lagrangian are

$L=\sum_{1\leq i<j\leq
d}\left|\Phi\partial_{i}\partial_{j}\Phi-\partial_{i}\Phi\partial_{j}\Phi\right|^{2}+V(|\Phi|^{2}).$
(64)

In the symmetry broken phase we would choose the potential in (64) to be a
Mexican hat. We can then expand the Lagrangian around one of the minima of the
potential, and to leading order in the fields we would arrive at the continuum
Lagrangian for the XY-plaquette model given in equation (18). This is a simple
check that the Lagrangian we wrote down has the desired symmetry broken phase
that we analysed in section 3.

To access the unbroken phase one would need to choose a generic potential with
a minimum at $|\Phi|^{2}=0$. Around this minimum the Lagrangian to lowest
orders reads

$L=\sum_{1\leq i<j\leq
d}\left|\Phi\partial_{i}\partial_{j}\Phi-\partial_{i}\Phi\partial_{j}\Phi\right|^{2}+m^{2}|\Phi|^{2}+\lambda|\Phi|^{4}.$
(65)

Unfortunately this effective Lagrangian is non-Gaussian even at lowest order
due to the kinetic term containing 4 $\Phi$ fields. Hence an analytic or even
perturbative description of the dynamics in the unbroken phase is beyond our
reach (as was also the case in Pretko_2018 .) This phase can still be studied
numerically, either on the lattice using (16) or in the continuum field theory
(65). This phase may also be accessible analytically using non-perturbative
techniques, like perhaps taking some sort of large $N$ limit.

It would also be of great interest to study the phase transition between
broken and unbroken phases. However, as the effective Lagrangian (65) is non-
Gaussian and strongly coupled in the vicinity of the phase transition, even a
mean field approximation of the transition is not possible. It is not even
clear if the phase transition is continuous or first order, and the literature
on phase transitions in similar models is uncertain.

For example the plaquette Ising model suzuki1972solution , which is the
$\mathbb{Z}^{2}$ version of the classical XY-plaquette model (16), has been
shown to undergo a weakly first order phase transition mueller2014nonstandard
; johnston2017plaquette . On the other hand a recent study of the phase
transition in models related to the $\mathbb{Z}^{N}$ versions of the X-cube
model indicated that the transition is continuous for $N>4$
lake2021subdimensional . As the plaquette Ising model is dual to the
$\mathbb{Z}^{2}$ X-cube model Vijay_2016 , these two results are in agreement.
One may then think of the XY-plaquette model as a $N\rightarrow\infty$ limit
of $\mathbb{Z}^{N}$ models, in which case the phase transition is likely
continuous, though more numerical or analytical evidence is needed to confirm
this conjecture.

###  5.2 Generalizations and interacting theories

Aside from modelling the unbroken phase and the phase transition, it would
also be interesting to construct interacting continuum field theories based on
the classical XY-plaquette model. We present two general systematic
constructions of such theories: the first method is to gauge the $U(1)$
subsystem symmetry to a tensor gauge theory in the spirit of Pretko_2018 ,
while the second idea is to promote a global $U(1)$ symmetry to a
subdimentional symmetry.

The most straightforward generalization of the classical XY-plaquette model is
to gauge the symmetry (19). This can be done by starting with the generalized
Lagrangian 64, and then replacing the derivative term
$\Phi\partial_{i}\partial_{j}\Phi-\partial_{i}\Phi\partial_{j}\Phi$ with a
“covariant” derivative term similar to Pretko_2018 ,

$D_{ij}[\Phi,\Phi]=\Phi\partial_{i}\partial_{j}\Phi-\partial_{i}\Phi\partial_{j}\Phi-
iA_{ij}\Phi^{2}.$ (66)

Here $A_{ij}$ is a symmetric $U(1)$ tensor gauge field with the transformation
properties $A_{ij}\rightarrow A_{ij}+\partial_{i}\partial_{j}\alpha$, while
$\Phi\rightarrow e^{i\alpha}\Phi$ under a gauge transformation.

We would also need to add a kinetic term for this gauge field, and in general
we one could investigate this tensor gauge theory in its own right. The gauge
invariant field strength of this theory would be Seiberg:2020wsg

$F_{[ij]k}=\partial_{i}A_{jk}-\partial_{j}A_{ik},$ (67)

and the kinetic term takes the form $F_{[ij]k}F^{[ij]k}$. We note that such a
theory only makes sense in $d\geq 3$, similar to how standard gauge theory
only makes sense in $d\geq 2$. From here one can continue studying
conservation laws, defects, and dynamics, similar to the gauge theories
considered in Seiberg:2020bhn ; Seiberg:2020wsg ; Pretko_2017 ; Pretko_2018 .

A different approach to systematically construct interacting theories from the
classical XY-plaquette model is by promoting a global $U(1)$ symmetry to a
subsytem symmetry. This procedure is reminiscent of the gauge principle, where
one promote a global symmetry to a local one.

For example we can start out with a field theory involving a single complex
scalar field $\psi$ which is charged under a global $U(1)$ symmetry, and a
Lagrangian

$L=\partial_{\mu}\psi^{\dagger}\partial^{\mu}\psi+V(\psi^{\dagger}\psi).$ (68)

Then we would like to partially gauge the global symmetry into a subsystem
symmetry. This can be done by introducing a “covariant” derivative

$D_{\mu}\psi=\left(\partial_{\mu}-ig\partial_{\mu}\phi\right)\psi,$ (69)

where $\phi$ acts as the analog of a gauge field. Under the subsystem symmetry
$\phi$ transforms as (19), while $\psi\rightarrow
e^{ig\sum_{\mu}f_{\mu}(x^{\mu})}\psi$. The full action of this theory would
then become

$L=D_{\mu}\psi^{\dagger}D^{\mu}\psi+V(\psi^{\dagger}\psi)+\sum_{\mu\neq\nu}\partial_{\mu}\partial_{\nu}\phi\partial^{\mu}\partial^{\nu}\phi.$
(70)

The resulting theory is an interacting field theory with a subsystem symmetry.

In essence one can partially gauge any theory with a conserved $U(1)$ charge,
resulting in a whole new family of theories to play with. It would be
interesting to investigate how these partially gauged theories behave, how to
generalize them to non-Abelian symmetries, and perhaps to understand if they
have analogous Coulomb and Higgs phases. Finally, we note that this partial
gauging procedure can also be used to construct interacting theories coupled
to the XY-plaquette models of Seiberg:2020bhn ; Seiberg:2020wsg , which would
also be interesting to investigate.

## Acknowledgments

We would like thank Michael Hermele, Nathan Seiberg, Shu-Heng Shao and Willy
Fischler for useful discussions. The work of AK and AR was supported, in part,
by the U.S. Department of Energy under Grant No. DE-SC0022021 and a grant from
the Simons Foundation (Grant 651440, AK). The work of JD was supported in part
by the National Science Foundation under Grant No. PHY–1914679.

##  Appendix A Correlation functions in the Classical XY-plaquette model

The structure of correlation functions of simple uncharged operators is very
non-standard in models with subsystem symmetry. In particular these
correlation functions tend to diverge even when the operators are far apart,
but one of the spatial variables is close gorantla2021lowenergy . Such UV
divergences in IR correlation functions has been considered a hallmark of
UV/IR mixing. As such, we are interested in computing the correlation
functions of the simple operators $\partial_{i}\partial_{j}\phi$, and to
observe any patterns of UV/IR mixing in the classical XY-plaquette model. We
will focus our analysis on the theory in two and three dimensions for
simplicity, and we think this also provides a general qualitative idea of the
structure of the correlation functions even in higher dimensions.

#### 2 dimensions:

For $d=2$ we can compute the correlation function directly as

$\left<\partial_{1}\partial_{2}\phi(x)\partial_{1}\partial_{2}(0)\phi\right>=g\int\frac{d^{2}k}{(2\pi)^{2}}e^{i(k_{1}x_{1}+k_{2}x_{2})}=g\delta(x_{1})\delta(x_{2}).$
(71)

This correlation function is somewhat simplistic, but does contain some
interesting features. Notice that the correlation function is indeterminate
when $x_{1}$ or $x_{2}$ are zero, which resembles the UV/IR mixing of
correlation functions in gorantla2021lowenergy . However, unlike in the XY
plaquette model of Seiberg:2020bhn ; gorantla2021lowenergy , this divergent
behavior cannot be regularized by an IR regulator (say by putting the system
on a torus,) but rather requires a UV regulator. A similar picture appears in
3 dimensions, though there the correlation functions have more structure.

#### 3 dimensions:

In three dimension we can compute the 2-point correlation function of
$\partial_{1}\partial_{2}\phi$ to be

$\displaystyle\left<\partial_{1}\partial_{2}\phi(x)\partial_{1}\partial_{2}(0)\phi\right>$
$\displaystyle=g\int\frac{d^{3}k}{(2\pi)^{3}}\frac{k_{1}^{2}k_{2}^{2}}{k_{1}^{2}k_{2}^{2}+k_{1}^{2}k_{3}^{2}+k_{2}^{2}k_{3}^{2}}e^{i(k_{1}x_{1}+k_{2}x_{2}+k_{3}x_{3})}$
(72)
$\displaystyle=\frac{g}{2}\int\frac{d^{2}k}{(2\pi)^{2}}\frac{|k_{1}||k_{2}|}{\sqrt{k_{1}^{2}+k_{2}^{2}}}e^{i(k_{1}x_{1}+k_{2}x_{2})-\frac{|k_{1}||k_{2}||x_{3}|}{\sqrt{k_{1}^{2}+k_{2}^{2}}}}$
$\displaystyle=\frac{g}{8\pi^{2}}\int d\theta
drr^{2}|\sin\theta\cos\theta|e^{ir(\cos\theta x_{1}+\sin\theta
x_{2})-r|\sin\theta\cos\theta||x_{3}|}$
$\displaystyle=\frac{g}{4\pi^{2}}\int_{0}^{2\pi}d\theta\frac{x_{3}\sin^{2}\theta\cos^{2}\theta\left(x_{3}^{2}\sin^{2}\theta\cos^{2}\theta-\left(x_{1}\cos\theta+x_{2}\sin\theta\right)^{2}\right)}{\left(x_{3}^{2}\sin^{2}\theta\cos^{2}\theta+(x_{1}\cos\theta+x_{2}\sin\theta)^{2}\right)^{3}}$

This result is finite if $x_{i}\neq 0$ for $i=1,2,3$, but diverges when the
two operators share one of the same coordinates. We can expand this integral
around these limits to see how it diverges. If we take $x_{3}\rightarrow 0$
then we see that

$\left<\partial_{1}\partial_{2}\phi(x_{1},x_{2},x_{3}\rightarrow
0)\partial_{1}\partial_{2}(0)\phi\right>=\frac{g\sqrt{x_{1}^{2}+x_{2}^{2}}}{8\pi|x_{1}x_{2}|}\cdot
x_{3}^{-2}+O(1).$ (73)

On the other hand if we take $x_{1}\rightarrow 0$ then

$\left<\partial_{1}\partial_{2}\phi(x_{1}\rightarrow
0,x_{2},x_{3})\partial_{1}\partial_{2}(0)\phi\right>=\frac{g}{8\pi
x_{3}^{2}}\cdot|x_{1}|^{-1}+O(1),$ (74)

and taking $x_{2}\rightarrow 0$ results in a similar expression with $x_{2}$
replacing $x_{1}$.

We can also work in cylindrical coordinates, defining
$x_{1}=\rho\cos(\delta)$, $x_{2}=\rho\sin(\delta)$, in which case the 2-point
function becomes

$\left<\partial_{1}\partial_{2}\phi(x)\partial_{1}\partial_{2}(0)\phi\right>=`\frac{g}{4\pi^{2}x_{3}^{3}}\int_{0}^{2\pi}d\theta\frac{\sin^{2}\theta\cos^{2}\theta\left(\sin^{2}\theta\cos^{2}\theta-\frac{\rho^{2}}{x_{3}^{2}}\cos^{2}(\theta-\delta)\right)}{\left(\sin^{2}\theta\cos^{2}\theta+\frac{\rho^{2}}{x_{3}^{2}}\cos^{2}(\theta-\delta)\right)^{3}},$
(75)

where the integral is only a function of $\rho/x_{3}$. In the limit
$x_{3}\ll\rho$ we get

$\left<\partial_{1}\partial_{2}\phi(x)\partial_{1}\partial_{2}(0)\phi\right>\approx\frac{g}{4\pi|\sin(2\delta)|x_{3}^{2}\rho},$
(76)

while in the limit $\rho\ll x_{3}$ the integral evaluates to

$\left<\partial_{1}\partial_{2}\phi(x)\partial_{1}\partial_{2}(0)\phi\right>\approx\frac{g|\sin\delta|+|\cos\delta|}{4\pi|\sin(2\delta)|x_{3}^{2}\rho},$
(77)

In all these limits the 2-point function diverges even though the operators
are separated. This is similar to the divergence of correlation functions as
two operators are brought close together, only in this case the operators stay
separated and only the separation in one of their coordinates is taken to be
close. This behavior is similar to the one obsercved in gorantla2021lowenergy
.

However, unlike in Seiberg:2020bhn ; gorantla2021lowenergy , this divergence
must be regulated by a UV regulator, and cannot be regulated by an IR
regulator (i.e. putting the system in a finite volume.) To see this, we can
take space to be a 3-Torus with equal lengths $\ell$, so that the momenta take
the discrete values $k_{i}=2\pi n_{i}/\ell$. Now the 2-point function becomes

$\displaystyle\left<\partial_{1}\partial_{2}\phi(x)\partial_{1}\partial_{2}(0)\phi\right>$
$\displaystyle=\frac{g}{\ell^{3}}\sum_{n_{1},n_{2},n_{3}\in\mathbb{Z}}\frac{n_{1}^{2}n_{2}^{2}}{n_{1}^{2}n_{2}^{2}+n_{1}^{2}n_{3}^{2}+n_{2}^{2}n_{3}^{2}}e^{i2\pi(n_{1}x_{1}+n_{2}x_{2}+n_{3}x_{3})/\ell}$
(78)

Taking $x_{1}=0$ it is clear that the sum over $n_{1}$ at fixed $n_{2},n_{3}$
does not converge, as is the case when taking $x_{2}=0$. If we take $x_{3}=0$
then we can still sum over $n_{3}$, giving us

$\displaystyle\left<\partial_{1}\partial_{2}\phi(x_{1},x_{2},0)\partial_{1}\partial_{2}(0)\phi\right>$
$\displaystyle=\frac{g}{\ell^{3}}\sum_{n_{1},n_{2},n_{3}\in\mathbb{Z}}\frac{n_{1}^{2}n_{2}^{2}}{n_{1}^{2}n_{2}^{2}+n_{1}^{2}n_{3}^{2}+n_{2}^{2}n_{3}^{2}}e^{i2\pi(n_{1}x_{1}+n_{2}x_{2})/\ell}$
(79) $\displaystyle=\frac{\pi
g}{\ell^{3}}\sum_{n_{1},n_{2}\in\mathbb{Z}}\frac{|n_{1}n_{2}|}{\sqrt{n_{1}^{2}+n_{2}^{2}}}e^{i2\pi(n_{1}x_{1}+n_{2}x_{2})/\ell}\coth\left(\frac{\pi|n_{1}||n_{2}|}{\sqrt{n_{1}^{2}+n_{2}^{2}}}\right),$

but now the sum over $n_{1}$ and $n_{2}$ is not convergent, and it still
requires a UV regulator.

As this divergence cannot be regulated by IR regulator, it is unclear to us
how these correlation functions fit into the picture of UV/IR mixing.
Furthermore, one can define a re-normalized operator
$\mathcal{O}_{12}\sim\Lambda^{3/2}\partial_{x}\partial_{y}\phi$, where
$\Lambda$ is some UV scale. Then this renormalized operator would have a
finite 2-point function of the form

$\left<\mathcal{O}_{12}(x_{1},x_{2},x_{3})\mathcal{O}_{12}(y_{1},y_{2},y_{3})\right>=\left\\{\begin{array}[]{cc}C,&x_{3}-y_{3}=x_{1}-y_{1}=0,\text{
or }x_{3}-y_{3}=x_{2}-y_{2}=0,\\\ 0,&\text{otherwise,}\end{array}\right.$ (80)

where $C$ is some constant that depends on the re-normalization scheme. This
renormalized operator has restricted mobility as it can only correlate along
lines of fixed $x_{1}$ and $x_{2}$ in the $x_{1}$–$x_{2}$ plane, and in terms
of these renormalized operators the theory would have no knowledge of the UV
physics.

However we are unconvinced that this is the correct way to describe the IR
theory, and are skeptical if this divergence of the 2-point function should be
regulated at all. Rather, we think this divergence is contains important
information about the IR physics, and the causal structure of the theory. An
analogous situation arises in relativistic quantum field theory, where it is
common for propagators diverge on the light cone. This light cone singularity
is a genuine property of the IR physics and the structure of space-time, and
not something which is needs to be regulated away. In any case it would be
interesting to further understand the consequences of these divergences, and
to fit them more squarely into the picture of UV/IR mixing.

## References

  * (1) A. Paramekanti, L. Balents and M. P. Fisher, _Ring exchange, the exciton bose liquid, and bosonization in two dimensions_ , _Physical Review B_ 66 (2002) 054526.
  * (2) C. Chamon, _Quantum Glassiness_ , _Phys. Rev. Lett._ 94 (2005) 040402 [cond-mat/0404182].
  * (3) J. Haah, _Local stabilizer codes in three dimensions without string logical operators_ , _Physical Review A_ 83 (2011) 042330.
  * (4) R. M. Nandkishore and M. Hermele, _Fractons_ , _Annual Review of Condensed Matter Physics_ 10 (2019) 295–313.
  * (5) M. Pretko, X. Chen and Y. You, _Fracton phases of matter_ , _International Journal of Modern Physics A_ 35 (2020) 2030003.
  * (6) N. Seiberg, _Field theories with a vector global symmetry_ , _SciPost Physics_ 8 (2020) .
  * (7) N. Seiberg and S.-H. Shao, _Exotic Symmetries, Duality, and Fractons in 2+1-Dimensional Quantum Field Theory_ , 2003.10466.
  * (8) N. Seiberg and S.-H. Shao, _Exotic $U(1)$ Symmetries, Duality, and Fractons in 3+1-Dimensional Quantum Field Theory_, _SciPost Phys._ 9 (2020) 046 [2004.00015].
  * (9) Y. You, J. Bibo, T. L. Hughes and F. Pollmann, _Fractonic critical point proximate to a higher-order topological insulator: How does uv blend with ir?_ , 2101.01724.
  * (10) Y. You and R. Moessner, _Fractonic plaquette-dimer liquid beyond renormalization_ , _arXiv preprint arXiv:2106.07664_ (2021) .
  * (11) P. Gorantla, H. T. Lam, N. Seiberg and S.-H. Shao, _The low-energy limit of some exotic lattice theories and uv/ir mixing_ , 2021.
  * (12) C. D. Batista and Z. Nussinov, _Generalized Elitzur’s theorem and dimensional reduction_ , _Phys. Rev. B_ 72 (2005) 045137 [cond-mat/0410599].
  * (13) S. Vijay, J. Haah and L. Fu, _Fracton topological order, generalized lattice gauge theory, and duality_ , _Physical Review B_ 94 (2016) .
  * (14) Y. You, T. Devakul, F. J. Burnell and S. L. Sondhi, _Subsystem symmetry protected topological order_ , _Phys. Rev. B_ 98 (2018) 035112.
  * (15) Y. You, Z. Bi and M. Pretko, _Emergent fractons and algebraic quantum liquid from plaquette melting transitions_ , _Physical Review Research_ 2 (2020) .
  * (16) M. Qi, L. Radzihovsky and M. Hermele, _Fracton phases via exotic higher-form symmetry-breaking_ , _Annals of Physics_ 424 (2021) 168360.
  * (17) E. Lake and M. Hermele, _Subdimensional criticality: condensation of lineons and planons in the x-cube model_ , 2021.
  * (18) D. Bulmash and M. Barkeshli, _Higgs mechanism in higher-rank symmetric u(1) gauge theories_ , _Physical Review B_ 97 (2018) .
  * (19) F. Strocchi, _Symmetry breaking_ , vol. 643. Springer, 2005.
  * (20) D. Gaiotto, A. Kapustin, N. Seiberg and B. Willett, _Generalized global symmetries_ , _Journal of High Energy Physics_ 2015 (2015) .
  * (21) E. Lake, _Higher-form symmetries and spontaneous symmetry breaking_ , _arXiv preprint arXiv:1802.07747_ (2018) .
  * (22) X.-C. Wu, W. Ji and C. Xu, _Categorical symmetries at criticality_ , _Journal of Statistical Mechanics: Theory and Experiment_ 2021 (2021) 073101.
  * (23) P. Gorantla, H. T. Lam, N. Seiberg and S.-H. Shao, _More Exotic Field Theories in 3+1 Dimensions_ , _SciPost Phys._ 9 (2020) 073 [2007.04904].
  * (24) A. Karch and A. Raz, _Reduced conformal symmetry_ , _Journal of High Energy Physics_ 2021 (2021) .
  * (25) M. Pretko, _The fracton gauge principle_ , _Physical Review B_ 98 (2018) .
  * (26) A. Gromov, _Towards classification of fracton phases: The multipole algebra_ , _Phys. Rev. X_ 9 (2019) 031035.
  * (27) M. Suzuki, _Solution and critical behavior of some” three-dimensional” ising models with a four-spin interaction_ , _Physical Review Letters_ 28 (1972) 507.
  * (28) M. Mueller, W. Janke and D. A. Johnston, _Nonstandard finite-size scaling at first-order phase transitions_ , _Physical Review Letters_ 112 (2014) 200601.
  * (29) D. A. Johnston, M. Mueller and W. Janke, _Plaquette ising models, degeneracy and scaling_ , _The European Physical Journal Special Topics_ 226 (2017) 749.
  * (30) M. Pretko, _Subdimensional particle structure of higher rank u(1) spin liquids_ , _Physical Review B_ 95 (2017) .