<figure><img src="image_1.png" /></figure>

HAL
open science

# Hollow lattice tensor gauge theories with bosonic matter

José M Cruz, Masafumi Udagawa, Pedro Bicudo, Pedro Ribeiro, Paul A Mcclarty

► To cite this version:

José M Cruz, Masafumi Udagawa, Pedro Bicudo, Pedro Ribeiro, Paul A Mcclarty. Hollow lattice tensor gauge theories with bosonic matter. 2025. hal-05406619

HAL Id: hal-05406619

https://hal.science/hal-05406619v1

Preprint submitted on 9 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Hollow Lattice Tensor Gauge Theories with Bosonic Matter

José M. Cruz,¹,²,³ Masafumi Udagawa,² Pedro Bicudo,¹ Pedro Ribeiro,¹,⁴ and Paul A. McClarty⁵,⁶

¹CeFEMA-LaPMET, Departamento de Física, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal

²Department of Physics, Gakushuin University, Mejiro, Tokyo, 171-8588, Japan

³Department of Physics, University College Dublin, Dublin, D04V1W8, Ireland*

⁴Beijing Computational Science Research Center, Beijing 100193, China

⁵Laboratoire Léon Brillouin, CEA, CNRS, Université Paris-Saclay, CEA-Saclay, 91191 Gif-sur-Yvette, France

⁶Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str. 38, 01187 Dresden, Germany†

Higher rank gauge theories are generalizations of electromagnetism where, in addition to overall charge conservation, there is also conservation of higher rank multipoles such as the total dipole moment. In this work we study a four dimensional lattice tensor gauge theory coupled to bosonic matter which has second rank tensor electric and magnetic fields and charge conservation on individual planes. Starting from the Hamiltonian, we derive the lattice action for the gauge fields coupled to $q = 1, 2$ charged scalars. We use the action formulation to carry out Monte Carlo simulations to map the phase diagram as a function of the gauge ($\beta$) and matter ($\kappa$) couplings. We compute the nature of correlators at strong and weak coupling in the pure gauge theory and compare the results to numerical simulations. Simulations show that the naive weak coupling regime (small $\kappa$, large $\beta$) does not survive in the thermodynamic limit. Instead, the strong coupling confined phase, spans the whole phase diagram. It is a proliferation of instantons that destroys the weak coupling phase and we show, via a duality transformation, that the expected strong confinement is present in the analog of Wilson line correlators. For finite matter coupling at $q = 1$ we find a single thermodynamic phase albeit with a first order phase transition terminating in a critical endpoint. For $q = 2$ it is known that the X-cube model with $\mathbb{Z}_2$ fractonic topological order is recovered deep in the Higgs regime. The simulations indeed reveal a distinct Higgs phase in this case.

## CONTENTS

<table><tr><td>I. Introduction</td><td>2</td></tr><tr><td>II. Hollow Rank-2 Gauge Theory with Bosonic Matter</td><td>2</td></tr><tr><td>A. Hamiltonian and Lattice Action</td><td>2</td></tr><tr><td>B. Remarks on the Monte Carlo method</td><td>4</td></tr><tr><td>C. Overview</td><td>4</td></tr><tr><td>III. Limiting Cases</td><td>5</td></tr><tr><td>A. Switching off the gauge dynamics</td><td>5</td></tr><tr><td>B. Infinite matter coupling</td><td>5</td></tr></table>

<table><tr><td>C. Infinite gauge coupling</td><td>6</td></tr><tr><td>D. Pure gauge theory at strong coupling</td><td>6</td></tr><tr><td>IV. Pure Gauge Theory: Duality, Instantons and Confinement</td><td>8</td></tr><tr><td>V. Phase Diagram of the A tensor gauge theory</td><td>11</td></tr><tr><td>A. Correlation functions and pinch points</td><td>15</td></tr><tr><td>VI. Discussion</td><td>17</td></tr><tr><td>A. Derivation of the lattice action</td><td>18</td></tr><tr><td>B. Pinch Point Correlations from the Continuum Field Theory</td><td>19</td></tr></table>

<footnote>* jose.m.cruz@ucdconnect.ie</footnote>

<footnote>† paul.mc-clarty@cnrs.fr</footnote>
C. Symmetries of Hollow Rank-2 U(1) 21

1. Correlation Functions 22

References 22

# I. INTRODUCTION

Condensed matter systems have long furnished physicists with a tremendous variety of distinct states of matter and it has been the work of many decades to find organizing principles general enough to encompass them. Following the work of Wilson [1] and Polyakov [2] on strongly coupled gauge theories on the lattice, it became necessary to generalize older ideas of symmetry breaking and local order parameters that were the staple of traditional solid state physics. Phases characterized by nonlocal order parameters and higher form symmetries arise naturally among certain classes of condensed matter systems such as quantum spin liquids where the low energy physics is captured by some kind of gauge principle.

Among such unusual states of matter are those with exotic-seeming conservation laws such as conservation of dipole moments [3, 4]. Consideration of models founded on such conservation laws came in the wake of the discovery of fracton phases [5–8] where quasiparticles arise with restricted mobility. For more details on fractons, together with references to the quite extensive literature, we refer the reader to reviews [9–12] on the subject. Following the discovery of such states of matter it was natural to look for gauge theories that are to these higher rank conservation laws as ordinary U(1) gauge theory is to charge conservation.

The resulting counterparts to ordinary Maxwell electrodynamics, higher rank gauge theories, have symmetric tensorial gauge fields and a generalized Abelian gauge invariance [3, 4, 13–24]. In addition to their roots in fracton phases they have connections to certain classical constrained systems with glassy dynamics and elastic media [11, 25, 26]. Their weak coupling continuum description also carries surprises by containing short wavelength low energy modes thus breaking the usual clean separation between UV and IR physics [20].

Our interest, in this paper, is to provide a global view of the phase diagram of lattice-regularized higher rank gauge theories in cases where the naive continuum limit breaks down [27]. To this end we provide a combined numerical and analytical study of the hollow rank-two theory in four dimensions [15, 28, 29], otherwise known as A tensor gauge theory [20] coupled to scalar matter. We present results of Monte Carlo simulations as a function of both gauge $\beta$ and matter $\kappa$ couplings for both charge one $q=1$ and charge two $q=2$ Higgs fields. The latter

case is of interest as a limit of the theory has fracton topological order [28, 29].

This work parallels early work on lattice Abelian Higgs models [30, 31] where the strong coupling pure gauge theory was shown to exhibit a confined phase. In 3+1-dimensions, in the pure U(1) gauge theory, there is a phase transition to a deconfined phase where Maxwell electromagnetism describes the low energy physics. The confined and deconfined phases of the pure gauge theory can be characterized by, respectively, an area and perimeter law for Wilson line operators. Coupling to scalar matter replaces the area law by a perimeter law because the Wilson line is screened by the charged matter. But the deconfined phase has power law correlations while the confined phase that smoothly connects to a Higgs regime is gapped.

For the pure gauge theory, we show that the finite size simulation results are compatible with weak and strong coupling expansions. However, finite size scaling reveals that features from the naive continuum limit at weak coupling do not survive the thermodynamic limit. As we know from earlier work [27], proliferation of instantons destroy the analog of the deconfined phase. We revisit this issue from a modern perspective by passing to dual variables. As we show in detail, the instantons that disorder the gauge fields lead to behavior that qualitatively matches the strong coupling behavior indicating a single phase for all $\beta$. This is the higher rank analog of the Polyakov confinement mechanism present in pure U(1) gauge theory in 2+1 dimensions and in two-form gauge theory in 3+1 dimensions [32]. Pinch points, that would be expected in certain electric field correlation functions from the continuum tensor gauge field theory [33], are absent in the quantum phase diagram.

From these considerations combined with simulation results we conclude that the Higgsed A tensor gauge theory with charge $q=1$ has a single phase passing from Higgs to confined regimes but with a line of first order phase transitions terminating at a critical endpoint reminiscent of a liquid-gas phase diagram.

In contrast, in the $q=2$ theory, there are two phases separated by a first order line. One of these phases – that we may call the Higgs phase – is continuously connected to the X-cube model [28, 29] and therefore has fractonic topological order.

# II. HOLLOW RANK-2 GAUGE THEORY WITH BOSONIC MATTER

## A. Hamiltonian and Lattice Action

A natural starting point to formulate the lattice gauge theory is first to introduce global symmetries associated with currents $J_0$ and $J_{ij}$, that is symmetric in its indices
i, j ∈ {x, y, z} and i ≠ j, satisfying continuity equations

$$ \partial_0 J_0 = \sum_{i<j} \partial_i \partial_j J_{ij} \qquad (1) $$

implying conserved charges on planes. We may gauge these global symmetries by introducing gauge fields $A_0 \to A_0 + \partial_0 \Omega$ and $A_{ij} \to A_{ij} + \partial_i \partial_j \Omega$. The minimal coupling in the action $A_0 J_0 + A_{ij} J_{ij}$ gauges these global symmetries.

We then pass to the lattice by placing $U(1)$ rotor variables on plaquettes of a three dimensional cubic lattice indexed by vertices $\vec{r}$. Spatial plaquettes or faces are indexed by $\vec{r}$ and $ij = xy, xz, yz$. On each spatial plaquette there is a phase variable $A_{ij}(\vec{r})$ and conjugate number operator denoted $E_{ij}(\vec{r})$ obeying

$$ [A_{ij}(\vec{r}), E_{ij}(\vec{r})] = i\delta_{\vec{r},\vec{r}'} (\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}). \qquad (2) $$

The lattice Gauss law is

$$ \sum_{i<j} \Delta_i \Delta_j E_{ij} = \rho. \qquad (3) $$

where $\rho$ is the scalar matter density defined on vertices and about which we say more below.

We further introduce operators $U_{ij}(\vec{r}) \equiv \exp(ia^2 A_{ij}(\vec{r}))$. Gauge transformations act on these as

$$ U_{xy}(\vec{r}) \to U_{xy}(\vec{r}) \eta(\vec{r}) \eta^{-1}(\vec{r} + \vec{e}_x) \qquad (4) \\ \eta(\vec{r} + \vec{e}_x + \vec{e}_y) \eta^{-1}(\vec{r} + \vec{e}_y) $$

where $\eta(\vec{r})$ is a phase shift acting on a vertex. This is the lattice version of the transformation of continuum $A_{xy}$ given above and which has a straightforward generalization to the other plaquette operators.

From these fields we may build gauge invariant local field strengths $L_{[ij]k}(\vec{r})$ with $i \neq j \neq k$:

$$ L_{[ij]k} = \exp(i(\Delta_i A_{jk} - \Delta_j A_{ik})) \qquad (5) $$

supported on elementary tubes on the lattice (as illustrated in the right panel of Fig. 1).

<figure><img src="image_2.png" /><figcaption>FIG. 1: Representation of operators $L_{[ij]k}$ and $L_{ij\tau}$ introduced in the main text.</figcaption></figure>

The Hamiltonian for the pure lattice gauge theory is

$$ H_{\text{Gauge}} = \sum_{\vec{r}, i<j} \frac{1}{g_e^2} E_{ij}^2 + \frac{1}{g_m^2} (L_{[ij]k}(\vec{r}) + L_{[ij]k}^{\dagger}(\vec{r})). \qquad (6) $$

In the literature, this theory goes by several names: the A tensor gauge theory [20], the rank-2 hollow $U(1)$ gauge theory [28], and the (0, 1)-scalar rank-2 gauge theory [29].

In addition we introduce scalar matter fields $\theta(\vec{r})$ coupled to the gauge fields together with conjugate fields $L(\vec{r})$:

$$ [\theta(\vec{r}), L(\vec{r}')] = i\delta_{\vec{r},\vec{r}''}. \qquad (7) $$

We choose the Higgs fields to have integer charge $q$ such that the Hamiltonian is

$$ H_{\text{Higgs}} = \sum_{\vec{r}} \frac{L(\vec{r})^2}{2M} - V \sum_{i<j} \cos(\Delta_i \Delta_j \theta(\vec{r}) + qA_{ij}(\vec{r})) \qquad (8) $$

and the Gauss law is, more precisely,

$$ \sum_{i<j} \Delta_i \Delta_j E_{ij} = qL(\vec{r}). \qquad (9) $$

The matter-gauge coupling term in the Hamiltonian is gauge invariant under simultaneous local shifts of the phase of $\theta(\vec{r})$ and on the corners of the $A_{ij}(\vec{r})$ plaquette operator. Evidently, application of the second term in Eq. 8 to the vacuum results in the creation of a plaquette excitation and alternating charges $\pm q$ around the plaquette. Successive application of such operators on neighbouring plaquettes leads to hopping of pairs of charges, forming a dipole, accompanied by a flux string.

The overall Hamiltonian is $H = H_{\text{Gauge}} + H_{\text{Higgs}}$. To gain insight into this model and for the purposes of investigating the strongly interacting regime numerically it will be convenient to write down a lattice action.

$$ S_{\text{Gauge}} = -\beta \sum_{\tau=1}^{N} \sum_{\vec{r}} \left( \cos(\Delta_i A_{jk}(\vec{r}) - \Delta_j A_{ki}(\vec{r})) + \cos(\Delta_i \Delta_j A_0^{\tau}(\vec{r}) - \Delta_0 A_{ij}^{\tau}(\vec{r})) \right) \qquad (10) $$

where the fields are defined on a hypercubic lattice indexed by $\vec{r} = (\hat{r}, \hat{x}, \hat{y}, \hat{z})$. The gauge invariant objects in this action are illustrated in Fig. 1. In addition to the gauge invariant space cubes $L_{[ij]k}$ (Fig. 1(left)) defined on a spatial cube, there are elementary time cubes $L_{ij\tau}$ (Fig. 1(right)) with pairs of spatial plaquettes connected by four temporal links and spatial cubes.

The Higgs-gauge part of the action is

$$ S_{\text{Higgs}} = -\kappa \sum_{\tau=1}^{N} \sum_{\vec{r}} \left( \cos(qA_0^{\tau}(\vec{r}) + \Delta_0 \theta^{\tau}(\vec{r})) + \sum_{i<j} \cos(qA_{ij}^{\tau}(\vec{r}) + \Delta_i \Delta_j \theta^{\tau}(\vec{r})) \right) \qquad (11) $$

The gauge invariant objects in the gauge-matter action are illustrated in Fig. 2. These are defined on temporal
links (left) and on spatial plaquettes (right) with matter
fields appearing at the vertices. Under gauge transforma-
tions, the action is left invariant by imposing $\theta \rightarrow \theta - q\Omega$.

<figure><img src="image_3.png" /></figure>

FIG. 2: Representation of gauge invariant combinations of gauge and matter fields. On the left we have the time link with matter fields at the ends $N_{\tau}$. On the right we have the space plaquette with matter fields on the corners $N_{ij}$.

A derivation of the action $S = S_{\text{Gauge}} + S_{\text{Higgs}}$ from the Hamiltonian is given in appendix A.

Somewhat similarly to the generalization of 1-form
to higher form gauge theories on the lattice, the rank-
two theory considered here can be generalized to higher
rank. For example, instead of plaquette (2-surface) terms
$A_{ij}(\vec{r})$ in the rank-two theory, the rank-three theory has
variables defined on 3-surfaces.

B. Remarks on the Monte Carlo method

In order to determine the nature of the phases in the
theory defined in the previous section, we performed
direct Monte Carlo simulations on the lattice action
(Eqs. 10 and 11). The Monte Carlo simulations were
performed using the Metropolis-Hastings algorithm on a
hypercubic $N^4$ lattice with periodic boundary conditions.
Efficient updating of the configuration was achieved by
using OpenMP and a division of the variables to be up-
dated into disjoint sets labelled by a parity. For tempo-
ral links this is a 2-parity, while for spatial plaquettes a
4-parity was used. For each such variable, the lattice is
divided into staples — the collection of variables on which
the update of a single plaquette, link or vertex variable
depends. We show the staples for the update of temporal
links (in Fig. 3), spatial plaquettes (Fig. 4) and matter
vertices (Fig. 5).

<figure><img src="image_4.png" /></figure>

FIG. 3: Time Staples - staples associated with the tem-
poral link variables. One such variable - the one to be
updated - is shown in yellow. Each time link is shared
between 12 time cubes $L_{ij\tau}$ (4 for each $ij$ plane) and by
one matter $N_{\tau}$ object.

<figure><img src="image_5.png" /></figure>

FIG. 4: Space Staples - staples associated with the spa-
tial plaquettes gauge variables. One such variable - the
one to be updated - is shown in yellow. Each spatial pla-
quette is shared between 4 space cubes $L_{[ij]k}$, one matter
$N_{ij}$ object and by 2 time cubes $L_{ij\tau}$.

<figure><img src="image_6.png" /></figure>

FIG. 5: Matter Staples - staples associated with the matter fields variables. One such variable - the one to be updated - is shown in yellow. Each matter field is shared between 12 $N_{ij}$ objects (4 for each $ij$ plane) and by 2 $N_{\tau}$ objects.

The autocorrelation time was monitored across the
simulations. The observed upper bound on the thermal-
ization time was 1500 Monte Carlo sweeps for the pa-
rameters studied. In practice, 10⁴ sweeps were chalked
up as a thermalization period and dropped from the be-
ginning of each Monte Carlo time series. The autocorre-
lation times in the thermalized time series were at most
200 sweeps for the parameters studied and measurements
were taken at intervals of 10³ sweeps to ensure that they
were sufficiently decorrelated. For error estimation all
error bars were computed using the jackknife method.

C. Overview

Let us step back and summarize what we have seen
so far. As emphasized in the introduction, we have
a lattice model with an Abelian gauge invariance of a
non-standard variety. Charges are conserved but iso-
lated charges are not present in the theory. Instead the
simplest charge sector consists of four charges that can
be moved around in pairs respecting the conservation
of charge in each plane of the lattice. The theory has
no Lorentz invariance in the continuum limit nor even
continuous rotational invariance and is instead tailored
to exotic condensed matter applications stemming from
fracton models such as the X-cube model [34].
The rest of this article is devoted to an exploration of the quantum phases obtained by tuning the matter couplings. In this context it will be illuminating to refer to the behavior of the ordinary lattice Abelian Higgs model in $(3+1)D$ with lattice action

$$S_{\text{ABH-Gauge}} = -\beta \sum_{\tau=1}^{N} \sum_{\vec{r}} \left( \begin{aligned} & \cos(\Delta_i A_j(\vec{r}) - \Delta_j A_i(\vec{r})) && (12) \\ & + \cos(\Delta_i A_0^\tau(\vec{r}) - \Delta_0 A_i^\tau(\vec{r})) \end{aligned} \right)$$

and

$$S_{\text{ABH-Higgs}} = -\kappa \sum_{\tau=1}^{N} \sum_{\vec{r}} \left( \begin{aligned} & \cos(qA_0^\tau(\vec{r}) + \Delta_0\theta^\tau(\vec{r})) && (13) \\ & + \sum_{i<j} \cos(qA_i^\tau(\vec{r}) + \Delta_i\theta^\tau(\vec{r})). \end{aligned} \right)$$

Here there is $A_i(\vec{r}) \to A_i(\vec{r}) + \Delta_i \Omega(\vec{r})$, $\theta(\vec{r}) \to \theta(\vec{r}) - \Omega(\vec{r})$ gauge invariance. While this model is well known we briefly review the salient features. In the pure gauge theory $\kappa = 0$, in $3+1$D the large $\beta$ phase is in the so-called Coulomb phase with a gapless, linearly dispersive photon and gapped magnetic monopole excitations [30, 35]. In this phase, electric charges are deconfined and the naive continuum field theory of scalar Maxwell electrodynamics describes the long wavelength behavior of the theory. In contrast, in the small $\beta$ phase the photon is gapped out. Separating electric charges produces an electric flux string between them and hence a linear confining potential. Relatedly Wilson line operators obey an area law in the confined phase and a perimeter law in the deconfined law. The Coulomb phase is stable to switching on finite $\kappa$ but there is a first order phase boundary at some $\kappa(\beta)$ at which the gapped charges condense into a Higgs phase and the photon is gapped out. From the pure gauge confined phase, screening by charges means that the area law is replaced by a perimeter law for infinitesimal coupling $\kappa$. Osterwalder and Seiler [36] and then Fradkin and Shenker [30] showed that there is a continuous region for $q=1$ where Green's functions are analytic that overlaps with parts of the Higgs phase and into the confining regime. These are therefore the same thermodynamic phase. When $q=2$ instead, the Higgs and confined regimes are different phases. The phase boundary between them separates a regime with deconfined $Z_2$ charges (Higgs) from a phase where all charges are gapped.

As we shall see the phase diagram for the A tensor gauge theory differs from its relative – the $3+1$ dimensional Abelian-Higgs model – in several important respects.

### III. LIMITING CASES

To gain an understanding of the nature of the phases of the lattice gauge theory it is useful, first, to consider

limiting cases where one coupling is either sent to zero or to infinity corresponding to the four bounding edges of the phase diagram.

#### A. Switching off the gauge dynamics

The simplest case of the four is the one where we switch off the gauge coupling $\beta \to 0$ entirely so that the action only has the $S_{\text{Higgs}}$ term. In the Hamiltonian language this is equivalent to sending $g_m/g_e \to \infty$. We may go to the unitary gauge $\theta^\tau(\vec{r}) = 0$ to reveal that the model consists only of decoupled plaquettes. The theory is therefore trivial in this limit.

#### B. Infinite matter coupling

If instead we take the matter-gauge coupling to infinity $\kappa \to \infty$, the matter fields and gauge fields are locked to one another with the constraint $pA_{ij}^\tau(\vec{r}) + \Delta_i \Delta_j \theta^\tau(\vec{r}) = 0$ with a similar expression for the temporal links. Once again it is easiest to see what happens in the unitary gauge. There one finds the constraint $A_{ij}^\tau(\vec{r}) = 2\pi n_{ij}^\tau(\vec{r})/q$ where $n_{ij}^\tau(\vec{r})$ is integer valued. In the following we consider two cases $q=1$ and $q=2$. In the former case, only gauge equivalent configurations appear in the theory for any $\beta$ so again the theory is trivial. However, in the $q=2$ case, the plaquette and link fields are $Z_2$ gauge fields and the theory is a pure fractonic $Z_2$ gauge theory.

We examine the $q=2$ case from the perspective of the Hamiltonian. The analog of the $\kappa \to \infty$ limit is then $MV \to \infty$. In this limit, the matter field ceases to be dynamical and we may gauge it away. With $V$ the only scale in the problem, the single plaquette ground states are those with $|A_{ij}(\vec{r}) = 0, \pi\rangle$ when $q=2$. Considering now the pure gauge Hamiltonian in this subspace we observe that the $L_{[ij]k}$ operators are diagonal in the single plaquette basis with eigenvalues $+1$ for $|A_{ij}(\vec{r}) = 0\rangle$ and $-1$ for $|A_{ij}(\vec{r}) = \pi\rangle$. In other words, we may replace these operators with products of $\sigma^z$ Pauli operators on the elementary spatial tubes. Raising and lowering operators between these states modulo $2\pi$ are $\exp(\pm i\pi E_{ij}(\vec{r}))$. Within the large $V$ subspace, these act like $\sigma^x$ on plaquettes. We may then introduce a term into the Hamiltonian that enforces the Gauss law dynamically. Since the matter fields have dropped out this takes the form

$$U \sum_{\vec{r}} \prod_{\{ij\} \in \vec{r}} (-)^{\Delta_i \Delta_j E_{ij}(\vec{r})} \sim U \sum_{\vec{r}} \prod_{\{ij\} \in \vec{r}} \sigma_{ij}^x(\vec{r}) \quad (14)$$

where the product runs over the 12 plaquettes connected to site $\vec{r}$. The dual cubic lattice has sites at the center of cubes on the direct lattice. Mapping the model to the dual lattice leads to variables living on links instead of plaquettes. The resulting model is recognizably the X-cube model which has fractonic excitations [34]. The observation that the X-cube model arises deep within the Higgs regime of the hollow rank-2 gauge theory coupled to $q=2$ scalar matter was reported in Refs. [28, 29].
In particular, the X-cube mode arises when the matter dynamics has been frozen out and where the gauge electric field term has been set to zero. When $g_e$ is finite, the gauge fields become dynamical. Within the 2D local subspace such dynamics can be generated by a term of the form

$$W \sum_{\vec{r},ij} \sigma_{ij}^{x}(\vec{r}). \qquad (15)$$

We expect the fracton phase, whose parent Hamiltonian is the X-cube model, to be stable to the presence of the dynamical term up to some threshold. Switching once again to the action, we therefore expect a phase transition from the fractonic phase at large $\beta$ to a paramagnetic phase at small $\beta$.

### C. Infinite gauge coupling

In the limit $\beta \to \infty$, the gauge dynamics is frozen out and we may gauge away the $A_{ij}(\vec{r})$ fields leaving the matter action

$$S_{\beta \to \infty} = - \sum_{\tau=1}^{N} \sum_{\vec{r}} \kappa \cos(\Delta_0 \theta^\tau(\vec{r})) + \sum_{i<j} \kappa \cos(\Delta_i \Delta_j \theta^\tau(\vec{r})). \qquad (16)$$

In the ordinary compact $U(1)$ Abelian Higgs model, this limit leads to a lattice $XY$ model in 4D with exchange on links and with critical coupling $\kappa_c$. In the rank two case there is an $XY$ degree of freedom on all sites but the temporal and spatial couplings are inequivalent. The latter is analogous to a ring exchange term and this difference results in a considerable enhancement of the symmetries. While the ordinary lattice $XY$ model has a single global $U(1)$ symmetry this $XY$-plaquette model has subsystem symmetries. To see this select one plane in a spatial direction (for example indexed by position $\vec{r}_0$ and directions $ij = xy, xz$ or $yz$) and perform a global shift on this plane of $\phi \to \phi + \alpha$. The action is invariant under this shift. There are $L_x + L_y + L_z$ such planes but there are constraints coming from the fact that a shift on all fixed $xy$ planes simultaneously is equivalent to a shift on all $xz$ or all $yz$ planes. In total there are $L_x + L_y + L_z - 2$ $U(1)$ subsystem symmetries. One expects spontaneous breaking of these symmetries at large $\kappa$ and a paramagnetic phase at small $\kappa$.

The continuum field theory that is expected to describe the long wavelength behavior of this model is

$$\mathcal{L} = \frac{m}{2} (\partial_0 \phi)^2 - \frac{1}{4\mu} (\partial_i \partial_j \phi)^2. \qquad (17)$$

We have carried out a Monte Carlo simulation of the model with action Eq. 16 finding the expectation values of both the temporal links terms and the plaquette terms. Results are shown for a $10^4$ lattice in Fig. 6. The simulation results provide evidence of a first order phase transition with $\kappa_c \approx 0.55$.

<figure><img src="image_7.png" /><figcaption>FIG. 6: Figures showing the temporal-matter (top) and spatial-matter (bottom) expectation values for the $XY$ plaquette model both showing a first order transition at $\kappa_c \approx 0.55$.</figcaption></figure>

### D. Pure gauge theory at strong coupling

The final limiting case, $\kappa = 0$, is also the most involved. This is the one where the matter-gauge coupling is taken to zero leaving only the pure gauge theory of Eq. 11. In this subsection, we explore the strong ($\beta \ll 1$) coupling limit leaving the naive weak coupling analysis to an Appendix ($\beta \gg 1$). This organization reflects a crucial feature of the phase diagram of this model namely that the naive weak coupling theory breaks down as a consequence of the proliferation of instantons.

We begin the strong coupling analysis with the partition function in the zero matter coupling limit,

$$Z = \int \prod_{\vec{x},\mu,\tau} dA_{\mu}(\vec{x}, \tau)e^{-S[A]} \qquad (18)$$

where $\mu \in \{12, 13, 23, 0\}$ and $S[A] = S_{\text{ABH-Gauge}} + S_{\text{ABH-Higgs}}$ is the action derived in the previous section. Let us consider arbitrary combinations of the elementary gauge invariant objects shown in Fig. 1 that we call schematically $\exp(iE_{ij}^\tau(\vec{r}))$ for the time-space cube and $\exp(iB_{[ij]k}^\tau(\vec{r}))$ for the spatial hollow cube. We wish to evaluate

$$\left\langle \prod_{E \in C_E} e^{iE} \prod_{B \in C_B} e^{iB} \right\rangle. \qquad (19)$$
To take the strong coupling limit, we may expand $exp(-\beta S)$ in powers of $\beta$. Inspecting Eq. 19 term by term in this expansion, one sees that the integrals over the gauge fields vanish unless the observables to be computed cancel with the conjugate term coming from products of

$$ \frac{1}{2} \left[ e^{i(\Delta_i A_{jk}(\vec{r}) - \Delta_j A_{ki}(\vec{r}))} + e^{-i(\Delta_i A_{jk}(\vec{r}) - \Delta_j A_{ki}(\vec{r}))} \right] \quad (20) $$

or

$$ \frac{1}{2} \left[ e^{i(\Delta_i \Delta_j A_0^\tau(\vec{r}) - \Delta_0 A_{ij}^\tau(\vec{r}))} + e^{-i(\Delta_i \Delta_j A_0^\tau(\vec{r}) - \Delta_0 A_{ij}^\tau(\vec{r}))} \right] \quad (21) $$

in the action. If set $C_E$ contains $N_E$ elements and $C_B$ contains $N_B$ elements we find that

$$ \left\langle \prod_{E \in C_E} e^{iE} \prod_{B \in C_B} e^{iB} \right\rangle = \exp \left[ (N_E + N_B) \ln \left( \frac{\beta}{2} \right) \right]. \quad (22) $$

It is instructive to consider some special cases of this result. One case is the expectation values of elementary plaquettes. At strong coupling the general result reveals that

$$ \langle L_{ij\tau} \rangle_{\text{SC}} = \langle L_{[ij]k} \rangle_{\text{SC}} = \frac{\beta}{2}. \quad (23) $$

A second useful observable is obtained by taking the product of time-space cubes along the temporal direction. This is the hollow rank-2 analog of the Polyakov loop and we call it the time tube $L_{ij}$:

$$ L_{ij}(\vec{r}) = \prod_{n=0}^{N_t-1} L_{ij\tau}(\vec{r}+n\vec{e}_{\tau}) = \exp(-i \sum_{n}^{N_t-1} \Delta_i \Delta_j A^{\tau}(\vec{r}+n\vec{e}_{\tau})) \quad (24) $$

One can think of this as a world line of two pairs of opposite charges separated by one link and in the strong coupling limit we find

$$ \langle L_{ij}(\vec{r}) \rangle_{\text{SC}} = \exp \left( N_{\tau} \ln \left( \frac{\beta}{2} \right) \right). \quad (25) $$

Finally we generalize this result to the case where the charge world lines are placed on the corners of a spatial plaquette of length $R_i$ in one direction and $R_j$ in another direction measured in lattice spacings. This gives a handle on the potential involved in separating the charges.

$$ \langle L_{ij}(\vec{r}; R_i, R_j) \rangle_{\text{SC}} = \exp \left( N_{\tau} R_i R_j \ln \left( \frac{\beta}{2} \right) \right). \quad (26) $$

This shows that there is a linear confining potential between pairs of dipoles and a quadratic potential for separating pairs of dipoles into the constituent charges.

Even more generally, consider a large hollow spatial cube with edge length $L$. Then the expectation of this in the small $\beta$ phase goes like $exp(-L^3 \ln(\beta/2))$. Similarly for a time tube with both temporal and spatial extents being equal to $L$.

To probe the potential between charges it is necessary to extend the time-cubes in the temporal direction, the time tube. But this object decays exponentially with the temporal span of the system, making it difficult to probe accurately using numerical methods. To verify the validity of strong coupling expansion it is useful the look at a simpler object, the extended time-cube $L_{ij\tau}^r(\vec{r})$, which is analogous to the extended time-tube but without the temporal elongation.

$$ \langle L_{ij\tau}^r \rangle_{\text{SC}} = \exp \left( r^2 \ln \left( \frac{\beta}{2} \right) \right). \quad (27) $$

<figure><img src="image_8.png" /><figcaption>FIG. 7: Normalized logarithm of the extended time cubes for a system size $V = 24x24x24$. Next to the data we plot $r^2$ for comparison.</figcaption></figure>

It is worth comparing these results to the behavior of ordinary $U(1)$ lattice gauge theory in 4D at strong coupling. As is well known, the theory at strong coupling lives in a confined phase with an area law Wilson loop and a linear potential between pairs of charges. In the rank-2 theory, rather than pairs of charges pairs of dipoles are confined in this limit and the analog of the area law Wilson line comes from a loop that separates pairs of dipoles in space. The Polyakov loops are quantitatively identical in the rank-2 and ordinary $U(1)$ theories. However, the form of the gauge invariant objects in the action indicates the natural analogs of Wilson loops in the higher rank theory occupy a 3-volume. At strong coupling these then obey a volume law.

<figure><img src="image_9.png" /><figcaption>FIG. 8: Figure showing expanded time tubes with spatial extent $r'$ used to measure the quadrupolar potential.</figcaption></figure>

In Appendix B, we discuss the naive weak coupling limit where the pure gauge theory on the lattice maps
to a tensor gauge field theory. We summarize the main results here and compare them to the strong coupling case. For example, we find at weak coupling

$$ \langle L_{ij\tau} \rangle_{\text{WC}} = \exp(-C/\beta) \quad (28) $$

to be contrasted with Eq. 23 at strong coupling. The constant C is a non-universal constant.

For $L_{ij}$ wrapping around in the temporal direction with spatial extent $r$, the weak coupling result is

$$ \langle L_{ij}(\vec{r}; r) \rangle_{\text{WC}} = \exp \left( - \frac{r N_t C'}{\beta} \right) \quad (29) $$

for some constant $C'$ and the strong coupling case ($\kappa = 0$) is

$$ \langle L_{ij}(\vec{r}; r) \rangle_{\text{SC}} = \exp \left( r^2 N_t \log \left( \frac{\beta}{2} \right) \right). \quad (30) $$

We note that switching on an infinitesimal $\kappa$ around the large $\beta$ phase leads to

$$ \langle L_{ij}(\vec{r}; r) \rangle_{\text{SC}, \kappa \neq 0} = \exp(-\zeta(2r^2 + 4N_t)). \quad (31) $$

This is analogous to the effect of small $\kappa$ on the strong coupling phase of ordinary U(1) gauge theory where an area law for Wilson lines ($\kappa = 0$) is replaced by a perimeter law ($\kappa \neq 0$).

As in the standard U(1) gauge theory, the rank two theory has a softer potential for the separation of pairs of dipoles in the weak coupling limit than in the strong coupling case where the $r^2$ dependence is analogous to area law behavior. As the distance scaling is qualitatively different in the two limits they cannot be interpolated by varying $\beta$. We would therefore expect a phase transition to separate the two regimes if they were stable phases. As we show in the next section, the strong coupling regime survives to large $\beta$ as a result of non-perturbative physics.

## IV. PURE GAUGE THEORY: DUALITY, INSTANTONS AND CONFINEMENT

A classic result in theoretical physics is the Polyakov confinement mechanism in three dimensional pure U(1) lattice gauge theory where instantons proliferate in the weakly coupled gauge theory and disorder the gauge field leading to a mass gap. In the analog of this physics in four dimensional pure gauge theory, one finds that magnetic monopole worldlines proliferate only below some finite coupling and there is a finite $\beta$ phase transition.

It turns out that the hollow rank 2 gauge theory in four dimensions behaves similarly to the three dimensional U(1) gauge theory. Although proliferation of instantons in this theory was observed to occur in the very early literature on this topic [27], we have revisited the dual description of the model to show explicitly that the instanton lead to confinement at large $\beta$ in the sense that strong coupling and dual weak coupling calculations coincide. The absence of a phase transition is then naturally

to be expected and we explore this issue numerically in the next section.

The perspective in this section connects to the work of Shao and Seiberg that discusses variants of continuum higher rank field theories and their duals [20, 21, 24]. In particular, our work refers to the lattice version of the $\tilde{A}$ and $\tilde{\phi}$ theories discussed in Ref [20].

We take a cubic lattice in four dimensions with sites $\vec{r}$ and directions labelled $\tau, x, y, z$. We introduce elementary fields $A_{ij}(\vec{r})$ defined to be symmetric in their indices $i, j = x, y, z$ but with $i \neq j$ and living on plaquettes and $A_{\tau}(\vec{r})$ living on temporal links. These are gauge fields transforming as

$$ A_{\tau} \rightarrow A_{\tau} + \Delta_{\tau} \alpha \quad (32) $$

$$ A_{ij} \equiv A_{ij} + \Delta_i \Delta_j \alpha. \quad (33) $$

We further define gauge invariant tensor electric and magnetic fields

$$ E_{ij} \equiv \Delta_{\tau} A_{ij} - \Delta_{i} \Delta_{j} A_{\tau} \quad (34) $$

$$ B_{[ij]k} \equiv \Delta_i A_{jk} - \Delta_j A_{ik} \quad (35) $$

where $E_{ij}$ is symmetric in its indices while $B_{[ij]k}$ is antisymmetric in the $i, j$ indices and $k \neq i, j$. These definitions imply that

$$ B_{[xy]z} + B_{[yz]x} + B_{[zx]y} = 0. \quad (36) $$

The partition function for the theory is

$$ Z = \int DA \exp \left( -\beta \sum_{r,i<j} \left( \cos(B_{[ij]k}(\vec{r})) + \cos(E_{ij}(\vec{r})) \right) \right). \quad (37) $$

Here the $A$ fields are defined from $-\pi$ to $\pi$. Without worrying about the range of $A$ we note that the electric and magnetic fields satisfy the analog of a Bianchi identity

$$ \Delta_{\tau} B_{[ij]k} = \Delta_{i} E_{jk} - \Delta_{j} E_{ik}. \quad (38) $$

The behavior of the theory in the large $\beta$ limit is not as in the naive continuum limit. Instead one must consider carefully the effects of the compactness of the gauge fields. To make these effects more apparent we first “Vilainize” the action. The general idea is to make use of the identity:

$$ \exp(-\beta(1 - \cos\theta)) = \sum_{m=-\infty}^{\infty} e^{im\theta} e^{-\beta I_m(\beta)} \quad (39) $$

where the modified Bessel function approximates to $e^{\beta e^{-m^2/4\beta}/\sqrt{2\pi\beta}}$ at large $\beta$ or weak coupling. Note that this explores the regime where there is a naive continuum limit of small $\theta$ fluctuations in the action.

The analogs of integers $m$ for the hollow rank two theory are $n_{[ij]k}$ for $B_{[ij]k}$ and $n_{\tau ij}$ for $E_{ij}$. Use of the above identity in the action leads to the interpretation of the
gauge potential as Lagrange multipliers that we integrate out to get constraint equations. These are

$$ \Delta_{\tau} n_{\tau yz} + \Delta_{x} (n_{[xy]z} + n_{[xz]y}) = 0 \quad (40) $$

$$ \Delta_{\tau} n_{\tau xz} + \Delta_{y} (n_{[yz]x} + n_{[yx]z}) = 0 \quad (41) $$

$$ \Delta_{\tau} n_{\tau xy} + \Delta_{z} (n_{[zx]y} + n_{[zy]x}) = 0 \quad (42) $$

and

$$ \Delta_i \Delta_j n_{\tau ij} = 0. \quad (43) $$

We shall now introduce integer fields $\phi_{[ij]k}$ defined on the centers of spatial cubes. These allow one to solve the constraint equations as follows:

$$ n_{[ij]k} \equiv -\Delta_{\tau} \phi_{[ij]k} \quad (44) $$

$$ n_{\tau ij} \equiv \Delta_k (\phi_{[ki]j} + \phi_{[kj]i}) . \quad (45) $$

Another way of proceeding is to naively write the cosine as $(\theta - 2\pi n)^2$ with a sum over $n$. Then we write the $n$ fields using the solution above plus integer versions of the electric and magnetic fields with an integer potential $h$. We may absorb the $h$ fields into the $A$ fields rendering them non-compact. In other words we define

$$ \tilde{E}_{ij} \equiv \Delta_{\tau} A_{ij} - \Delta_{i} \Delta_{j} A_{\tau} + \Delta_{\tau} h_{ij} - \Delta_{i} \Delta_{j} h_{\tau} \quad (46) $$

$$ \tilde{B}_{[ij]k} \equiv \Delta_i A_{jk} - \Delta_j A_{ik} + \Delta_i h_{jk} - \Delta_j h_{ik} \quad (47) $$

and then we evaluate the cross terms in the square that appear between the $\tilde{E}$, $\tilde{B}$ and the $\phi$ terms. It turns out that, after partial integration, the terms vanish as a consequence of the Bianchi identity. So what about the remaining terms in $\phi$? Using either approach, these are:

$$ S_{\phi} = \frac{1}{4\beta} \sum_{r,i<j} \left[ (\Delta_{\tau} \phi_{[ij]k})^2 + (\Delta_k (\phi_{[ki]j} + \phi_{[kj]i}))^2 \right]. \quad (48) $$

The constraint on the components should be maintained when summing over configurations.

Expanding out the action and going to momentum gives the following kernel coupling $(\phi_{[xy]z}, \phi_{[yz]x}, \phi_{[zx]y})$:

$$ \mathbf{M}(\mathbf{k}) = \begin{pmatrix} \omega^2 - k_x^2 - k_y^2 & k_y^2 & k_x^2 & k_z^2 \\ k_y^2 & \omega^2 - k_y^2 - k_z^2 & k_z^2 & k_x^2 - k_z^2 \\ k_x^2 & k_z^2 & \omega^2 - k_x^2 - k_z^2 & k_z^2 \end{pmatrix} \quad (49) $$

The denominator in the inverse of this kernel is

$$ \omega^2 (\omega^4 - 2\omega^2(k_x^2 + k_y^2 + k_z^2) + 3(k_x^2 k_y^2 + k_x^2 k_z^2 + k_y^2 k_z^2)) \quad (50) $$

that gives the spectrum reported in the Seiberg et al. paper on their $\hat{\phi}$ theory [20].

Now we go back and add a general source term to the action of the form

$$ \sum_{\vec{r},i<j} J_{ij} A_{ij} + J_{\tau} A_{\tau} \quad (51) $$

so that now when the gauge fields are integrated out the constraint equations take the form

$$ \Delta_{\tau} n_{\tau yz} + \Delta_{x} (n_{[xy]z} + n_{[xz]y}) = J_{yz} \quad (52) $$

$$ \Delta_{\tau} n_{\tau xz} + \Delta_{y} (n_{[yz]x} + n_{[yx]z}) = J_{xz} \quad (53) $$

$$ \Delta_{\tau} n_{\tau xy} + \Delta_{z} (n_{[zx]y} + n_{[zy]x}) = J_{xy} \quad (54) $$

and

$$ \Delta_i \Delta_j n_{\tau ij} = J_{\tau}. \quad (55) $$

summed over $i < j$.

The analog of the continuity equation for this current is
$$ \Delta_{\tau} J_{\tau} = \Delta_{i} \Delta_{j} J_{ij} \quad (56) $$

again summed over $i < j$. As may be verified by looking at Eqs. 54.

In principle one may look for a general solution to these equations for the integer-valued $n$ and $\phi$. But we are interested in making contact with the strong coupling expansion. Therefore, we choose a particular form for the Wilson line mimicking the form we studied in that limit. This means choosing current lines $J_{\tau}(x, y, z, \tau) = \pm 1$ along four lines along the time direction centered at:

$$ J_{\tau} = +1 : x = -R/2, y = -R/2, z = 0 \quad (57) $$

$$ J_{\tau} = -1 : x = +R/2, y = -R/2, z = 0 \quad (58) $$

$$ J_{\tau} = -1 : x = -R/2, y = +R/2, z = 0 \quad (59) $$

$$ J_{\tau} = +1 : x = +R/2, y = +R/2, z = 0 \quad (60) $$

and zero everywhere else. This straightforwardly satisfies the continuity equation.

One may formally solve the equations as follows for this particular source term

$$ n_{\tau xy} \equiv \Delta_z (\phi_{[zx]y} + \phi_{[zy]x}) + (\Delta_x \Delta_y)^{-1} J_{\tau} \quad (61) $$

$$ n_{\tau xz} \equiv \Delta_y (\phi_{[yx]z} + \phi_{[yz]x}) \quad (62) $$

$$ n_{\tau yz} \equiv \Delta_x (\phi_{[xy]z} + \phi_{[xz]y}) \quad (63) $$

$$ n_{[ij]k} \equiv -\Delta_{\tau} \phi_{[ij]k} \quad \text{for } ij = xy, xz, yz \quad (64) $$

Now let us recall the Poisson summation formula that introduces continuous fields into the problem

$$ \sum_{p} f(p) = \sum_{q} \int d\psi f(\psi) e^{2\pi iq\psi} \quad (65) $$

Let us also neglect the constraint on the sums of the $\phi_{[ij]k}$.

We find, for the full dual action in the presence of currents:

$$ S_{\psi,a} = \frac{1}{4\beta} \sum_{r,i<j} (\Delta_{\tau} \psi_{[ij]k})^2 + (\Delta_k (\psi_{[ki]j} + \psi_{[kj]i}) + I_{ij\tau})^2 + 2\pi i a_{[ij]k} \psi_{[ij]k} \quad (66) $$

where

$$ I_{ij\tau} = \delta_{i,x} \delta_{j,y} (\Delta_x \Delta_y)^{-1} J_{\tau}. \quad (67) $$

We introduce the three component field $(\psi_{[xy]z}, \psi_{[yz]x}, \psi_{[zx]y})$ and defining the kernel

$$ \mathbf{M} \equiv \begin{pmatrix} \Delta_{\tau}^2 + \Delta_x^2 + \Delta_y^2 & -\Delta_y^2 & -\Delta_x^2 \\ -\Delta_y^2 & \Delta_{\tau}^2 + \Delta_y^2 + \Delta_z^2 & -\Delta_z^2 \\ -\Delta_x^2 & -\Delta_z^2 & \Delta_{\tau}^2 + \Delta_x^2 + \Delta_z^2 \end{pmatrix}. \quad (68) $$
It is important to note that the order of indices does not follow a cyclic permutation. The choice recorded above with $\psi_{[zx]y}$ instead of $\psi_{[xz]y}$ is consistent with the sign structure in the double derivative operator $M$.

Now we may integrate out $\psi$ leaving an action for the $a_{[ij]k}$ and the currents. We find

$$ \beta \left( 2\pi i a_{[ij]k} + \frac{1}{2\beta} \mathcal{I}_{[ij]k} \right) \mathbf{M}^{-1} \left( 2\pi i a_{[ij]k} + \frac{1}{2\beta} \mathcal{I}_{[ij]k} \right) \quad (69) $$

where, for the particular charge configuration of interest:

$$ (\mathcal{I}_{[xy]z}, \mathcal{I}_{[yz]x}, \mathcal{I}_{[zx]y}) = (0, \Delta_z (\Delta_x \Delta_y)^{-1} J_\tau, -\Delta_z (\Delta_x \Delta_y)^{-1} J_\tau) \quad (70) $$

To progress we regularize the potential writing

$$ \mathbf{M}^{-1}(\vec{r}, \vec{r}') = V_0 \delta_{\vec{r}, \vec{r}'} + V(\vec{r} - \vec{r}') \quad (71) $$

and we decouple the interaction between the $a_{[ij]k}$ fields by introducing new variables $\eta_{[ij]k}$.

We end up with the action

$$ S_{\eta,a} = \sum_{r,i<j} -\frac{1}{4\pi^2} \frac{1}{4\beta} \eta_{[ij]k} \mathbf{M} \eta_{[ij]k} + i a_{[ij]k} \eta_{[ij]k} \quad (72) \\ - 4\pi^2 \beta V_0 a_{[ij]k} a_{[ij]k} + 2i\pi a_{[ij]k} \mathbf{M}^{-1} \mathcal{I}_{[ij]k} $$

Let us examine the part of the action that depends on $a_{[ij]k}$. These are integer valued fields that we must sum over:

$$ \sum_{a_{[ij]k}=-\infty}^{\infty} e^{-4\pi^2 \beta V_0 a_{[ij]k} a_{[ij]k} + i a_{[ij]k} \eta_{[ij]k} + 2i\pi a_{[ij]k} \mathbf{M}^{-1} \mathcal{I}_{[ij]k}} \quad (73) $$

where clearly the regularization of the potential plays the role of suppressing the contribution from the $a_{[ij]k}$ when $\beta$ is sufficiently large (deep in the weak coupling regime). Essentially the monopole charges are then dilute and we may restrict the sum to $-1, 0, 1$. Doing this we find

$$ 1+2 \exp(-4\pi^2 \beta V_0) \cos(\eta_{[ij]k} + 2\pi \mathbf{M}^{-1} \mathcal{I}_{[ij]k}) + \dots \\ \approx \exp[2 \exp(-4\pi^2 \beta V_0) \cos(\eta_{[ij]k} + 2\pi \mathbf{M}^{-1} \mathcal{I}_{[ij]k})]. \quad (74) $$

Overall we find for the partition function

$$ \mathcal{Z}[\mathcal{L}] = \prod_{\vec{r}} \prod_{i<j \neq k} \int d\eta_{[ij]k} \exp \left[ \sum_{r,i<j} -\frac{1}{4\pi^2} \frac{1}{4\beta} \eta_{[ij]k} \mathbf{M} \eta_{[ij]k} \right. \\ \left. \qquad + 2 \exp(-4\pi^2 \beta V_0) \cos(\eta_{[ij]k} + 2\pi \mathbf{M}^{-1} \mathcal{I}_{[ij]k}) \right] \quad (75) $$

We may now compute the saddle point equations for the dynamical $\eta_{[ij]k}$ fields:

$$ \mathbf{M}_{ab}\eta_b = -M^2 \sin(\eta_a + 2\pi \mathbf{M}^{-1} \mathcal{I}_{[ij]k}) \quad (76) $$

where $a = [ij]k$. So there are three coupled differential equations here to solve in the presence of a source term.

We may solve these equations numerically without any further assumptions to understand the quantitative effective potential of the model at weak coupling. By introducing a fictional temperature $T$ we can use simulated annealing Monte Carlo to find a minimum for the action present in the exponent of 75.

Below we present the annealing results for a system with size $V = 30 \times 30 \times 30$ and the current distribution from equations 60.

<figure><img src="image_10.png" /><figcaption>FIG. 9: Right: Current configuration used. The red plane consists of positive charge and the blue plane of negative charge Center: $\eta_{[yz]x}$ field z-dependence for x and y inside the current configuration. The vertical lines correspond to the location of the charge planes of the corresponding color. Right: $\eta_{[yz]x}$ fields distribution for z=1</figcaption></figure>
The results from the annealing process show that the $\eta_{[yz]x}$ fields are accumulating near the currents distribution. As the distance from the places increases the field magnitude decreases exponentially, as can be seen from the fit in the center plot of figure 9. Outside of the planes the fields contributions are negligible. Similar results were obtained for the $\eta_{[zx]y}$, but with opposite sign. The contributions from $\eta_{[xy]z}$ are in general negligible when compared to the other fields.

The accumulation of the fields over the current distribution seems to suggest an area law. To verify this, the ground state energy obtained from the annealing procedure was calculated for different charge distributions with varying linear separations $R$ between charges.

<figure><img src="image_11.png" /><figcaption>FIG. 10: Ground state energy dependence on the charge pair separation $R$. The data was fitted to $a + bx + cx^2$ and the resulting curve was plotted on top of the data. The fitting parameters obtained were $a = 0.299$, $b = -5.02 \times 10^{-6}$ and $c = 1.18 \times 10^{-5}$, which shows the linear term is negligible and the expression obtained is purely quadratic in nature.</figcaption></figure>

By fitting the data we found that the ground state energy scales exactly like $R^2$, in agreement with the result found for the strong coupling regime.

This result implies that the point-like instanton solutions will not shield the strong potential from the strong-coupling regime and it will carry on for weak-coupling. Therefore there is no phase transition between these two regimes for the hollow rank-2 gauge theory.

## V. PHASE DIAGRAM OF THE A TENSOR GAUGE THEORY

We begin our exploration of the Higgsed A tensor lattice gauge theory by surveying the full phase diagrams of

the models specified by lattice action Eqs. 10 and 11 with charge $q = 1$ and 2 Higgs fields. In particular, we present results from Monte Carlo simulations as a function of $\beta$ and $\kappa$ computing the analog of the (gauge invariant) Polyakov loop operator

$$ \langle L_{ij} \rangle = \left\langle \prod_{n=0}^{N_t-1} L_{ij\tau} (\vec{r} + n\vec{e}_\tau) \right\rangle. \qquad (77) $$

Fig. 11 shows the variation of $\langle L_{ij} \rangle$ obtained from simulations on a $16^4$ hypercubic lattice for $q = 1$. Broadly, $\langle L_{ij} \rangle$ is strongly suppressed for small $\beta$ and finite for larger $\beta$. One finds that $\log\langle L_{ij} \rangle/N_t$ collapses results for different system sizes. This quantity is plotted along cuts (indicated as dashed lines in Fig. 11) in Fig. 12. Here one sees that there is an abrupt step indicative of a first order transition along cuts 2-5. The height of the step decreases as $\beta$ decreases and no such step is visible along cut 6 which, instead, is smooth as a function of the couplings. Therefore the existence of a critical point can be surmised. The white solid line in Fig. 11 marks the extent of the first order line and the approximate location of the putative critical endpoint.

<figure><img src="image_12.png" /><figcaption>FIG. 11: Phase diagram of the 4D charge 1 Higgsed hollow rank 2 gauge theory parametrized by $\beta$ and $\kappa$. The order parameter plotted is the time tube $\langle L_{ij} \rangle$. The small $\kappa$, large $\beta$ phase is connecting to the other phase by horizontal and vertical phase boundaries. There is a first order phase boundary (white line) connecting to a critical endpoint (white point). The color gradient is the time tube expectation value. Dashed black lines indicate lines along which cuts are shown in figures below. This phase diagram was obtained from Monte Carlo simulations on a $16^4$ lattice.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>(a) Line 1</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Line 2</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) Line 3</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) Line 4</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(e) Line 5</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(f) Line 6</figcaption></figure>

FIG. 12: Cuts represented in the time tubes phase diagram for charge $q = 1$ above. The inset in each image is the quotient of the time tubes logarithm with the system temporal span $N_τ$. This quantity is expected to be system size invariant. For strong coupling(small $β$) this sometimes fails (Lines 1, 4, 5). The expected value for the time tubes can be so small in this regime that it is unfeasible to obtain enough data for the necessary precision. For weak coupling(large $β$) it is always observed.
A cursory examination of Fig. 11 raises the possibility of a further transition, visible in $\langle L_{ij} \rangle$, for small fixed $\kappa$ as $\beta$ is varied. Looking instead at the horizontal cut 1 in Fig. 12 no sharp feature is visible. This requires further investigation including a careful finite size scaling analysis. We considered various observables finding no clear signs that a $\beta$ driven transition is present in the thermodynamic limit. In contrast, we find a compelling case for a crossover. For example, Fig. 13 shows data for the susceptibility of the logarithm of $L_{ij}$ plotted against $\beta/L^x$ for both zero matter coupling, $\kappa$ and small finite $\kappa = 0.5$ (the latter corresponding to cut 1 in Fig. 12). We find a collapse of data for different system sizes for exponent $x = 0.173$. This result implies that the peak in this quantity drifts to larger $\beta$ as the system size increases so the naive weak coupling theory does not survive in the thermodynamic limit.

<figure><img src="image_19.png" /><figcaption>FIG. 13: Susceptibility of the logarithm of the time tubes for $\kappa = 0$ (top) and $\kappa = 0.5$ (bottom). For both cases $x \approx 0.173$. A crossover between different regimes can be observed, but it does not survive the thermodynamic limit.</figcaption></figure>

The conclusion from this study is that there is a single

thermodynamic phase in the $q=1$ case but that there is a first order phase transition line similar to a liquid-gas transition line with a critical endpoint. This feature is also present in both the Abelian Higgs model in four dimensions (albeit accompanied by a deconfined phase for small $\kappa$, large $\beta$) and also in non-Abelian Higgs models (where, in common with the model considered here, there is a single phase). It would be of interest to investigate further the nature of the criticality at the endpoint but this is beyond the scope of the current investigation.

Now we turn to the charge $q=2$ case, where a similar study was employed. Fig. 14 shows the variation of $\langle L_{ij} \rangle$ obtained from simulations on a $16^4$ hypercubic lattice for $q=2$. We again can see that $\langle L_{ij} \rangle$ is strongly suppressed for small $\beta$ and finite for larger $\beta$. The main difference here is the first order transition extends until infinite $\kappa$ and lack of critical point. This transition can be clearly seen in the cuts 2-4 (indicated as dashed lines in Fig. 14) which are plotted in Fig. 15. Cuts 1 and 5 do not cross this transition and have a smooth crossover, similar to the $q=1$ case. We have already discussed the infinite matter coupling limit in section III B and concluded it should have a phase transition. We can then surmise that the observed first order transition extends until infinity. The white solid line in Fig. 14 serves as a guide for the eye for the first order transition location. The expected location of the X-Cube limit is also marked.

<figure><img src="image_20.png" /><figcaption>FIG. 14: Phase diagram of the 4D charge 2 higgsed hollow rank 2 gauge theory parametrized by $\beta$ and $\kappa$. The order parameter plotted is the time tube $\langle L_{ij} \rangle$. Unlike the charge 1 case, there is a phase boundary extending to infinite $\kappa$ separating the Higgs phase from the fractonic confined phase. This phase diagram was obtained from Monte Carlo simulations on a $16^4$ lattice.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>FIG. 15: Cuts represented in the time tubes phase diagram for charge $q = 2$ above. The inset in each image is the quotient of the time tubes logarithm with the system temporal span $N_τ$. This quantity is expected to be system size invariant. For strong coupling(small $β$) this sometimes fails for the same reason as before, but is always observed for weak coupling(large $β$)</figcaption></figure>
It can be concluded then that for $q = 2$ there is now an additional thermodynamic phase: a distinct Higgs phase analogous to that of the Abelian Higgs model in four dimensions, where a similar behavior is observed for charge $q = 2$ as well. However, the rank-2 Higgs phase has fractonic topological order.

As a first attempt to characterize this newfound fractonic Higgs phase we study the extended time cubes introduced in the previous section at different points in the phase diagram.

<figure><img src="image_22.png" /><figcaption>FIG. 16: Normalized logarithm of the extended time cubes for a system size $V=24 \times 24 \times 24$ in different locations of the phase diagram for charge $q = 1$. We have fitted the data to quadratic order and plotted the result. The fitting results were $0.034r^2 + 0.34r + 0.57$, $0.039r^2 + 0.046r + 0.99$ and $0.032r^2 + 0.005r + 0.97$ for the blue, orange and green curves respectively.</figcaption></figure>

As can be seen from the previous figure, there were no qualitatively different results for the extended timecubes in the fractonic Higgs phase, being well described by a quadratic function. This motivates the study of correlation functions as a way to characterize the different phases, which we do in the following section.

### A. Correlation functions and pinch points

So far we have found evidence for a single phase in the $q = 1$ hollow rank-two theory coupled to bosonic matter and for two phases in the $q = 2$ case. In this section, we consider the problem of characterizing these phases through correlation functions.

In particular, we compute the higher-rank electric tensor correlation function

$$C_{ij}^{kl}(r) = \langle L_{ij\tau}(r) L_{kl\tau}^{\dagger}(0) \rangle - \langle L_{ij\tau} \rangle \langle L_{kl\tau}^{\dagger} \rangle. \quad (78)$$

Analytically, these functions can be studied in the strong coupling regime similarly to the time tubes discussed earlier. To first-order in the strong coupling expansion there is no separation dependence, making it necessary to include higher order terms. This is done by finding the smallest possible gauge invariant combination that cancels every gauge variable and includes both the electric tensor we want to study. One example of such a construction can be found in Fig. 17.

<figure><img src="image_23.png" /><figcaption>FIG. 17: Strong coupling higher-order gauge variable structure for the correlation function $C_{ij}^{ij}$ when the electric tensor are separated in the $k$ direction. In the figure we represent the time links as dots on the corners of the space plaquettes. Since the electric field tensors are temporal-spatial objects, we present the two different time slices in which they live. The electric tensors in question are represented with a yellow dot.</figcaption></figure>

There is a $\beta/2$ factor for each individual time cube $L_{ij\tau}$ and space cube $L_{[ij]k}$ that is contained in the large structure that joins the electric tensors. For the case presented in Fig. 17 it is simple to see that the correlation function at distance $r$ to leading order in the strong coupling expansion is

$$C_{ij}^{ij}(r\vec{e}_k) = 2^{|r|} \left(\frac{\beta}{2}\right)^{4|r|} . \quad (79)$$

There are different possible combinations of $L_{ij\tau}$ and $L_{[ij]k}$ that lead to the same results. Those are taken into account with the factor $2^{|r|}$. Results for other correlation functions of interest are summarized in the table below:

<table><tr><td>$C_{ij}^{ij}(r\vec{e}_i/\vec{e}_j)$</td><td>$2^{|r|-1}(r-1)\left(\frac{\beta}{2}\right)^{4(|r|+2)}$</td></tr><tr><td>$C_{ik}^{kj}(r\vec{e}_k)$</td><td>$2^{|r|-1}(r-1)\left(\frac{\beta}{2}\right)^{4(|r|+2)}$</td></tr><tr><td>$C_{ik}^{kj}(r\vec{e}_i/\vec{e}_j)$</td><td>$2^{|r|}\left(\frac{\beta}{2}\right)^{4(|r|+1)}$</td></tr></table>

These correlation function were also obtained nu-
merically using the same Monte-Carlo algorithm used
throughout this work. Results are shown in Fig. 18 for
coupling $\beta = 1$ and $\kappa = 0$ and for $C_{xy}^{xy}$. Here we chose an
intermediate coupling as the absolute magnitude of the
correlation function is smaller for stronger couplings such

that it harder to distinguish from noise for $\beta$ much below
1. We observe a four-fold symmetry in the $xy$ plane and
two-fold symmetry in the $yz$ and $xz$ planes. The decay
along the $z$ direction is slow enough that we may fit the
result to the strong coupling prediction (Eq. 78) finding
good agreement for small distances.

<figure><img src="image_24.png" /><figcaption>FIG. 18: Logarithm of the correlation functions $C_{xy}^{xz}$ in real space for charge $q = 2$ and couplings $\beta = 1$ and $\kappa = 0$. The top plots are along the $z = 0$, $x = 0$ and $y = 0$ planes. The bottom cuts are along the $x$, $y$ and $z$ axis. For the latter we compare the numerical results with the strong coupling result and find agreement between them. For the former two this comparison is not made because the correlations decay too quickly (which is expected from the strong coupling expansion).</figcaption></figure>

Similar findings were found for the $C_{xy}^{xz}$ correlator including agreement with the strong coupling expansion.

Figs. 19 and 20 show the $C_{xy}^{xy}$ correlation function in momentum space for the $q = 2$ model at three different points in the phase diagram. The two plots on the left are for the confining phase respectively at intermediate and weak coupling. The plot on the right in each figure is taken within the Higgs phase. We have discussed the strong coupling result in real space. In momentum space the intensity is concentrated around the $p_x - p_y$ plane with a four-fold symmetry but with few features in that plane.

The case of weak coupling deserves further discussion.
The naive continuum limit of non-hollow rank two theory
reveals electric tensor correlation functions with pinch
point singularities [33] with four-fold symmetry. This
is reminiscent of the $U(1)$ spin liquid coming from the

deconfined phase of a pure $U(1)$ gauge theory which
also has pinch-points but with two-fold symmetry.
Pinch point provide a sharp experimental signature that
could be used to diagnose tensor Coulomb spin liquids in,
for example, neutron scattering experiments. We have
shown in Appendix B that similar features are present in
the naive weak coupling limit on the lattice in the theory
under consideration here.

We have also presented evidence that instanton prolif-
eration destroys the naive weak coupling regime. How-
ever, this evidence from the Monte Carlo required careful
finite size scaling. Inspection of the correlations at small
$\kappa$ and $\beta = 2$ that would exhibit pinch points in the naive
weak coupling regime instead show weak signs of a ridge
along the $p_x$ and $p_y$ axes which is suggestive of the pinch
point albeit only as a very weak feature visible here only
on a logarithmic scale. It is expected that this feature
scales away in the thermodynamic limit.
<figure><img src="image_25.png" /><figcaption>FIG. 19: Logarithm of the correlation functions $C_{xy}^{xy}$ in momentum space for charge $q = 2$ for different points in the phase diagram. In the left we have $\beta = 1, \kappa = 0$, in the middle $\beta = 2, \kappa = 0.5$ and on the right $\beta = 2, \kappa = 1.5$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>FIG. 20: Same plots as the figure above, but now with focus on the plane $p_z = 0$.</figcaption></figure>

The results presented above are expected to carry over to the phase of the $q = 1$ model. We now turn our attention to the other distinct phase: the Higgs phase in the $q = 2$ theory. Here the correlations again are concentrated in a plane but now the correlations in that plane show a distinctive “cross” feature along the $p_x$ and $p_y$ axes that is not present in the confined phase (see Fig. 20 — right panel).

## VI. DISCUSSION

Inspired by the continuum field theory counterpart that has been extensively studied over the last few years, in this paper we have carried out a systematic study of the hollow rank-two gauge theory in four dimensions on a lattice.

We have discussed the appropriate lattice Hamiltonian and action from first principles introducing generalizations of the plaquette variables one sees in one-form lattice gauge theory. We have also introduced scalar matter into the theory carrying charge $q = 1, 2$.

We have analysed various limiting cases of the hollow rank two Higgs theory. One interesting case is the limit of weak gauge coupling $\beta \to \infty$ controlled by matter-gauge coupling $\kappa$. We saw that the resulting statistical

mechanics model has $U(1)$ subsystem symmetries and, via Monte Carlo simulation, we observed that there is a first order transition at some $\kappa^*$ separating the paramagnet from a long-range ordered phase that spontaneously breaks subsystem symmetry.

We also studied the case where the matter coupling goes to zero. This pure gauge theory has a strong confinement regime at strong coupling that can be studied using perturbation theory where pairs of dipoles are confined by a linear potential. The usual continuum field theory for which we study the lattice analog has a number of peculiar properties that we describe in an appendix including pinch point correlations. However, this weak coupling picture is destroyed by instantons through a higher rank version of the Polyakov mechanism. To our knowledge, this work is the first to show, through a numerical solution of the Debye-Hückel equations, that strong confinement survives to weak coupling in this theory via a derivation of effective action for the instantons. We further bolster this result by studying the electric field correlations that show pinch points in the continuum weak coupling gauge theory. As expected no pinch points are present.

These findings for the limiting cases of the $q = 1$ theory reveal a single phase transition. Through Monte Carlo simulations of the full theory we find that the transition extends to finite $\beta, \kappa$ terminating at a critical endpoint that would be interesting to explore further. In other
words, this theory has a single phase that we characterized numerically and in the zero matter coupling setting.

The $q = 2$ theory has two phases. This can be seen already by taking the $\kappa \to \infty$ limit where, similarly to the charge two $U(1)$ gauge theory in four dimensions, there is a residual global $\mathbb{Z}_2$ degree of freedom. In this limit there is a $\beta$ driven phase transition and Monte Carlo simulations of the full theory show that this phase boundary meets the phase boundary originating from the infinite $\beta$ limit. This large-$\beta$, $\kappa$ Higgs phase is the dynamical gauge theory extension of the X-cube model.

There are various potentially interesting open problems that are beyond the scope of this study. It would be interesting to explore the full rank-2 scalar charge theory with diagonal terms as the mechanism of partial confinement that leads to the hollow theory studied here [28]. Also it would be interesting to see, via numerical simulation, that the weak coupling continuum theory survive in higher dimensions. More importantly, can confinement be avoided in four dimensions coupling to fermionic matter? What is the behavior on a lattice of other types of rank-$N$ theory? Recently, some of us revealed the presence of boundary phase transitions in generic Abelian higher form and non-Abelian gauge theories coupled to scalar matter [37]. What becomes of this transition in the case of higher rank gauge theories?

## ACKNOWLEDGEMENTS

PR and PB acknowledge support by FCT under contracts: LA/P/0095/2020, LaPMET, Laboratory of Physics for Materials and Emerging Technologies; and UID/04540, CeFEMA, Centre of Physics, Engineering and Advanced Materials. PM acknowledges funding from the CNRS. MU is supported by KAKENHI Grants No. JP22H01147 and No. JP23K22418 from JSPS of Japan.

## Appendix A: Derivation of the lattice action

In this appendix, we derive a lattice action for Hamiltonian $H = H_{\text{Gauge}} + H_{\text{Higgs}}$ given in the main text and repeated here for convenience.

$$H_{\text{Gauge}} = \sum_{\vec{r}, i<j} \frac{1}{g_e^2} E_{ij}^2 + \frac{1}{g_m^2} (L_{[ij]k}(\vec{r}) + L_{[ij]k}^{\dagger}(\vec{r})). \quad (A1)$$

and

$$H_{\text{Higgs}} = \sum_{\vec{r}} \frac{L(\vec{r})^2}{2M} - V \sum_{i<j} \cos(\Delta_i \Delta_j \theta(\vec{r}) + qA_{ij}(\vec{r})) \quad (A2)$$

with Gauss law constraint

$$G(\vec{r}) = \sum_{i<j} \Delta_i \Delta_j E_{ij}(\vec{r}) = qL(\vec{r}). \quad (A3)$$

Periodic boundary conditions are assumed, natural units will be used throughout and lattice constant $a$ is set to 1.

Here we implement the Gauss law by restrict to the constrained subspace in the derivation of the action through projector

$$P = \int D\phi \exp \left( i \sum_{\vec{r}} G(\vec{r}) \phi(\vec{r}) \right) \\ = \int D\phi \exp \left( i \sum_{\vec{r}} \sum_{i<j} \Delta_i^b \Delta_j^b E_{ij}(\vec{r}) \phi(\vec{r}) \right). \quad (A4)$$

where $D\phi = \prod_{\vec{r}} \frac{d\phi(\vec{r})}{2\pi}$

The partition function, in the physical subspace, is then

$$Z = \text{Tr} [\exp(-\beta H) P] = \text{Tr} \left[ \left( \exp \left( -\frac{\beta}{N} H \right) P \right)^N \right] \\ = \int DA \prod_{\tau=1}^{N} \langle A_{\tau} | \exp(-\Delta_{\tau} H) P | A_{\tau+1} \rangle \quad (A5)$$

where a time slicing has been carried out with $\Delta\tau = \beta/N$.

We then split the Hamiltonian into two parts $H = H_A + H_E$, where $H_A$ contains the terms that depend on the gauge variables $A_{ij}$, $H_E$ depending on the conjugate variables $E_{ij}$. Therefore, for small $\Delta\tau$ we have

$$\exp(-\Delta\tau H) = \exp(-\Delta\tau H_A - \Delta\tau H_E) \\ \approx \exp(-\Delta\tau H_A) \exp(-\Delta\tau H_E). \quad (A6)$$

Then we insert resolutions of the identity in the basis of $E$ eigenstates at each time slice

$$Z \approx \int DA \prod_{\tau=1}^{N} \langle A_{\tau} | \exp(-\Delta\tau H_A) \exp(-\Delta\tau H_E) P | A_{\tau+1} \rangle \\ = \int DA \prod_{\tau=1}^{N} \sum_{E_{\tau}} \langle A_{\tau} | \exp(-\Delta\tau H_A) | E_{\tau} \rangle \\ \qquad \langle E_{\tau} | \exp(-\Delta\tau H_E) P | A_{\tau+1} \rangle \\ = \int DAD\phi \sum_{E_{\tau}} \prod_{\tau=1}^{N} \exp(-\Delta\tau H_{A_{\tau}} - \Delta\tau H_{E_{\tau}}) \\ \qquad - i \sum_{\vec{r}} \sum_{i<j} \Delta_i^b \Delta_j^b E_{ij}^{\tau} \phi_{\tau}(\vec{r}) \\ \qquad + i \sum_{i<j, \vec{r}} A_{ij}^{\tau}(\vec{r}) E_{ij}^{\tau}(\vec{r}) - A_{ij}^{\tau+1}(\vec{r}) E_{ij}^{\tau}(\vec{r})$$
Rearranging and setting $\theta_{\tau}(\vec{r}) = -A_{0}^{\tau}(\vec{r})$ we find

$$Z = \text{const} \int DA \sum_{E_{\tau}} \exp\left(\sum_{\tau=1}^{N} \sum_{\vec{r},i<j} \frac{\Delta\tau}{g_m^2} \left(L_{[ij]k}^{\tau}(\vec{r}) + L_{[ij]k}^{\tau\dagger}(\vec{r})\right) \right. \\ \left. - \frac{\Delta\tau}{g_e^2} E_{ij}^{\tau}(\vec{r}) E_{ij}^{\tau\dagger}(\vec{r})\right. \\ \left. + i(A_0^{\tau}(\vec{r}) - A_0^{\tau}(\vec{r}+\mathbf{i}) - A_0^{\tau}(\vec{r}+\mathbf{j})) \right. \\ \left. + A_0^{\tau}(\vec{r}+\mathbf{i}+\mathbf{j}) + A_{ij}^{\tau}(\vec{r}) - A_{ij}^{\tau+1}(\vec{r})\right) E_{ij}^{\tau}(\vec{r}) \right)$$

where we have made a redefinition of the measure $DA_{\tau} \equiv DAD\phi$.

Now using

$$\exp(z \cos(\phi)) \approx \sum_n \exp\left(-\frac{1}{2z}n^2 + i\phi n\right). \quad (A7)$$

to carry out the sum over electric field variables we find

$$S = - \sum_{\tau, \vec{r}, i<j} \frac{2\Delta\tau}{g_m^2} \cos(A_{jk}(\vec{r} + \vec{i}) - A_{jk}(\vec{r})) \\ - A_{ki}(\vec{r} + \vec{j}) + A_{ki}(\vec{r}) \quad (A8) \\ + \frac{g_e^2}{2\Delta\tau} \cos(A_0^\tau(\vec{r}) - A_0^\tau(\vec{r} + \vec{i}) - A_0^\tau(\vec{r} + \vec{j})) \\ + A_0^\tau(\vec{r} + \vec{i} + \vec{j}) + A_{ij}^\tau(\vec{r}) - A_{ij}^{\tau+1}(\vec{r})$$

While this is the most general form of the action, we take the couplings constants for each therm to be identical $\frac{2\Delta\tau}{g_m^2} = \frac{g_e^2}{2\Delta\tau} \equiv \beta$ giving the gauge action in the main text. A similar calculation maybe carried out for the Higgs coupling leading to the final form of the lattice action.

## Appendix B: Pinch Point Correlations from the Continuum Field Theory

In the main text, we discuss at length the pure gauge theory in the strong coupling limit as this is the limit most relevant to characterising the gauge theory in 4D for all $\beta$. It is interesting, though, to contrast these results with those obtained by taking a naive continuum limit where the lattice plaquette length scale $a \to 0$.

Starting with lattice action:

$$S = -\beta \sum_{\tau} \sum_{\vec{r}} \left( \cos(\Delta_i A_{jk}(x) - \Delta_j A_{ki}(x)) \right. \\ \left. + \cos(\Delta_i \Delta_j A_0(x) - \Delta_0 A_{ij}(x)) \right). \quad (B1)$$

we obtain

$$S = \frac{\beta}{2} \sum_{\tau} \sum_{\vec{r}} (\Delta_i A_{jk}(x) - \Delta_j A_{ki}(x))^2 \\ + (\Delta_i \Delta_j A_0(x) - \Delta_0 A_{ij}(x))^2 \quad (B2)$$

to second order in the gauge fields.

We choose a gauge in which the “divergence-free” condition is satisfied at each vertex,

$$\begin{aligned} & - (A_0(x^\mu) - A_0(x^\mu - e^0)) + (A_{yz}(x^\mu) \\ & - A_{yz}(x^\mu - e^2) - A_{yz}(x^\mu - e^3) + A_{yz}(x^\mu - e^2 - e^3)) \\ & + (A_{zx}(x^\mu) - A_{zx}(x^\mu - e^3) \\ & - A_{zx}(x^\mu - e^1) + A_x(x^\mu - e^3 - e^1)) \\ & + (A_{xy}(x^\mu) - A_{xy}(x^\mu - e^1) \\ & - A_{xy}(x^\mu - e^2) + A_{xy}(x^\mu - e^1 - e^2))) = 0 \end{aligned} \quad (B3)$$

This is the direct generalization of Feynman gauge in the ordinary $U(1)$ case to this model.

Fourier transforming via

$$\left\{ \begin{aligned} A_{jk}(p_{\mu}) &= \frac{1}{\sqrt{N^4}} \sum_{p_{\mu}} e^{-ip_{\mu}x^{\mu}} A_{jk}(x^{\mu}) \\ A_0(p_{\mu}) &= \frac{1}{\sqrt{N^4}} \sum_{p_{\mu}} e^{-ip_{\mu}x^{\mu}} A_0(x^{\mu}) \end{aligned} \right. . \quad (B4)$$

the gauge constraint (B3) becomes

$$\begin{aligned} & - (1 - e^{-ip_0}) A_0(p_\mu) + (1 - e^{-ip_2})(1 - e^{-ip_3}) A_{23}(p_\mu) \\ & + (1 - e^{-ip_3})(1 - e^{-ip_1}) A_{31}(p_\mu) \\ & + (1 - e^{-ip_1})(1 - e^{-ip_2}) A_{12}(p_\mu) = 0. \end{aligned} \quad (B5)$$

The quadratic action is then

$$S = \frac{\beta}{2} [A_0(-p_{\mu}) \quad A_{23}(-p_{\mu}) \quad A_{31}(-p_{\mu}) \quad A_{12}(-p_{\mu})] \\ D^{-1}(p_{\mu}) \begin{bmatrix} A_0(p_{\mu}) \\ A_{23}(p_{\mu}) \\ A_{31}(p_{\mu}) \\ A_{12}(p_{\mu}) \end{bmatrix} \quad (B6)$$

with propagator

$$\begin{aligned} [D^{-1}(p_{\mu})]_{00} &= |e^{ip_1} - 1|^2 |e^{ip_2} - 1|^2 + |e^{ip_2} - 1|^2 |e^{ip_3} - 1|^2 \\ &\quad + |e^{ip_3} - 1|^2 |e^{ip_1} - 1|^2 \\ [D^{-1}(p_{\mu})]_{10} &= -(e^{ip_2} - 1)(e^{ip_3} - 1)(e^{-ip_0} - 1) \\ [D^{-1}(p_{\mu})]_{11} &= |e^{ip_0} - 1|^2 + 2|e^{ip_1} - 1|^2 \\ [D^{-1}(p_{\mu})]_{20} &= -(e^{ip_3} - 1)(e^{ip_1} - 1)(e^{-ip_0} - 1) \\ [D^{-1}(p_{\mu})]_{21} &= -(e^{-ip_2} - 1)(e^{ip_1} - 1) \\ [D^{-1}(p_{\mu})]_{22} &= |e^{ip_0} - 1|^2 + 2|e^{ip_2} - 1|^2 \\ [D^{-1}(p_{\mu})]_{30} &= -(e^{ip_1} - 1)(e^{ip_2} - 1)(e^{-ip_0} - 1) \\ [D^{-1}(p_{\mu})]_{31} &= -(e^{-ip_3} - 1)(e^{ip_1} - 1) \\ [D^{-1}(p_{\mu})]_{32} &= -(e^{-ip_3} - 1)(e^{ip_2} - 1) \\ [D^{-1}(p_{\mu})]_{33} &= |e^{ip_0} - 1|^2 + 2|e^{ip_3} - 1|^2 \end{aligned} \quad (B7)$$

We define its eigenvalues $\varepsilon_m(p_{\mu})$ and eigenvectors $\vec{u}^{(m)}(p_{\mu})$ by the usual relation,

$$D^{-1}(p_{\mu}) \vec{u}^{(m)}(p_{\mu}) = \varepsilon_{m}(p_{\mu}) \vec{u}^{(m)}(p_{\mu}). \quad (B8)$$
Out of the 4 eigenmodes, $m = 0, 1, 2, 3$, the lowest mode ($m = 0$) gives the zero mode: $\varepsilon_0(p_\mu) = 0$, as it corresponds to the motion of the gauge fields in the gauge-equivalent space. The corresponding eigenvector can be obtained from Eq. (B5) as

$$ \vec{u}^{(0)}(p_{\mu}) = \begin{bmatrix} A_0(p_{\mu}) \\ A_{23}(p_{\mu}) \\ A_{31}(p_{\mu}) \\ A_{12}(p_{\mu}) \end{bmatrix} = \begin{bmatrix} -(1 - e^{ip_0}) \\ (1 - e^{ip_2})(1 - e^{ip_3}) \\ (1 - e^{ip_3})(1 - e^{ip_1}) \\ (1 - e^{ip_1})(1 - e^{ip_2}) \end{bmatrix} . \quad (B9) $$

Expressions for the other eigenvectors cannot be found explicitly. We write the normal modes,

$$ \phi_m(p_\mu) = u_0^{(m)*}(p_\mu)A_0(p_\mu) + u_{23}^{(m)*}(p_\mu)A_{23}(p_\mu) \\ + u_{31}^{(m)*}(p_\mu)A_{31}(p_\mu) + u_{12}^{(m)*}(p_\mu)A_{12}(p_\mu), \quad (B10) $$

among which $\phi_0(p_\mu)$ vanishes identically, due to the gauge condition, Eq. (B5). Then the gauge fields can be expanded in the physical space as

$$ \begin{cases} A_0(p_\mu) = \sum_{m=1}^{3} u_0^{(m)}(p_\mu) \phi_m(p_\mu) \\ A_{ij}(p_\mu) = \sum_{m=1}^{3} u_{ij}^{(m)}(p_\mu) \phi_m(p_\mu) \end{cases} . \quad (B11) $$

With the normal modes, we can write down the action as

$$ S = \frac{\beta}{2} \sum_{p_{\mu}} \sum_{m=1}^{3} \varepsilon_m |\phi_m(p_{\mu})|^2 = \beta \sum_{p_{\mu}>0} \sum_{m=1}^{3} \varepsilon_m |\phi_m(p_{\mu})|^2. \quad (B12) $$

We are now in a position to compute correlation functions of operators. We focus on the time space cube

$$ \begin{aligned} L_{ij\tau}(x^{\mu}) &\equiv U_{\tau}(x^{\mu}) U_{\tau}^{-1}(x^{\mu} + e^{i}) U_{\tau}^{-1}(x^{\mu} + e^{j}) U_{\tau} \\ &\quad (x^{\mu} + e^{i} + e^{j}) U_{\tau}(x^{\mu}) U_{\tau}^{-1}(x^{\mu} + e^{0}) \quad (B13) \\ &= \exp \left[ i \left\{ \begin{aligned} &A_0(x^{\mu}) - A_0(x^{\mu} + e^{i}) - A_0(x^{\mu} + e^{j}) \\ &+ A_0(x^{\mu} + e^{i} + e^{j}) + (A_{23}(x^{\mu}) - A_{23}(x^{\mu} + e^{0})) \end{aligned} \right\} \right] \end{aligned} $$

and, for concreteness, we take $(i, j) = (2, 3)$ (the results can be generalized for other $(i, j)$ by symmetry arguments). In momentum space this is

$$ L_{23\tau}(x^{\mu}) = \exp \left[ i \frac{1}{\sqrt{N^4}} \sum_{p_{\mu}} \sum_{m=1}^{3} e^{ip_{\mu}x^{\mu}} \left\{ (1 - e^{ip_2})(1 - e^{ip_3}) u_0^{(m)}(p_{\mu}) \phi_m(p_{\mu}) + (1 - e^{ip_0}) u_{23}^{(m)}(p_{\mu}) \phi_m(p_{\mu}) \right\} \right]. \quad (B14) $$

We find:

$$ \begin{aligned} \langle L_{23\tau}(x^{\mu}) \rangle &= \exp \left[ -\frac{1}{N^4} \sum_{p_{\mu}, m} \frac{|(1 - e^{ip_2})(1 - e^{ip_3}) u_0^{(m)}(p_{\mu}) + (1 - e^{ip_0}) u_{23}^{(m)}(p_{\mu})|^2}{\beta \varepsilon_m(p_{\mu})} \right] \\ &\equiv \exp(-B/\beta). \quad (B15) \end{aligned} $$

and

$$ \begin{aligned} \langle L_{23\tau}^{\dagger}(x^{\mu}) L_{23\tau}(0) \rangle = & \exp \left[ -\frac{1}{N^4} \sum_{p_{\mu}} \sum_{m=1,2,3} \frac{1}{2\beta\varepsilon_m} |1 - e^{ip_{\mu}x^{\mu}}|^2 \right. \\ & \left. \qquad |(1 - e^{ip_2})(1 - e^{ip_3}) u_0^{(m)}(p_{\mu}) + (1 - e^{ip_0}) u_{23}^{(m)}(p_{\mu})|^2 \right]. \end{aligned} $$

It follows that the time tube expectation value goes like

$$ \langle L_{23}(x^{\mu}) \rangle = \exp(-N_{\tau}B/\beta). \quad (B16) $$

and the extended time tube goes like

$$ \langle L_{23}(x^{\mu}; r) \rangle = \exp \left( - \frac{rN_{\tau}B}{\beta} \right) \quad (B17) $$

We now consider the connected correlator

$$ C_{23\tau}(x^{\mu}) = \langle L_{23\tau}^{\dagger}(x^{\mu}) L_{23\tau}(0) \rangle - \langle L_{23\tau}(0) \rangle^2. \quad (B18) $$

which, in momentum space, is

$$ C_{23\tau}(p^{\mu}) = \frac{1}{\beta} \sum_{m=1,2,3} \frac{1}{\varepsilon_m} |(1 - e^{ip_2})(1 - e^{ip_3}) u_0^{(m)}(p_{\mu}) + (1 - e^{ip_0}) u_{23}^{(m)}(p_{\mu})|^2 \quad (B19) $$

We plot $C_{23\tau}(p^{\mu})$ for different momentum slices.

Fig. 21 shows the $p_2, p_3$ plane for fixed $p_0 = p_1 = 0$. In this plane there is a four-fold symmetry with a singular line along the $p_2 = 0$ and $p_3 = 0$ lines. Shifting the plane by plotting instead at fixed $p_1 = \pi/16$ reveals that where there were singular features there are valleys along these same lines (Fig. 22). The pinch points are localized in this plane.

<figure><img src="image_27.png" /><figcaption>FIG. 21: Four-fold pinch points for $p_0 = p_1 = 0$ and plotted in the $p_2, p_3$ plane.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>FIG. 22: Plane parallel to the four-fold pinch points for $p_0 = 0$, $p_1 = \pi/16$ and plotted in the $p_2, p_3$ plane showing how the pinch points are broadened into valleys.</figcaption></figure>

## Appendix C: Symmetries of Hollow Rank-2 U(1)

We now verify the symmetries associated with the hollow rank-2 U(1) action

$$ S_{\text{Gauge}} = -\beta \sum_{\tau=1}^{N} \sum_{\vec{r}} \left( \cos(\Delta_i A_{jk}(\vec{r}) - \Delta_j A_{ki}(\vec{r})) + \cos(\Delta_i \Delta_j A_0^\tau(\vec{r}) - \Delta_0 A_{ij}^\tau(\vec{r})) \right). \quad (C1) $$

First let us verify how it behaves with rotations. The figure below sketches how do the plaquette elements of the action vary with 90° rotations along some axis. The “orientation” of each plaquette is represented with an arrow.

<figure><img src="image_29.png" /><figcaption>FIG. 23: Changes due to 90° rotation around the k axis</figcaption></figure>

From the figure it is clear that rotations around the k axis affect plaquettes in the following manner

$$ A_{ij} \rightarrow A_{ij}, \quad A_{ik} \rightarrow -A_{jk}, \quad A_{jk} \rightarrow A_{ik}, \quad (C2) $$

where we need to take into the account the cyclic nature

of the indices when replacing them by actual coordinates $ijk = \{xyz, yzx, zxy\}$.

The relative positions between the colored dots also make it clear that the spatial differences need to transform like

$$ \Delta_k \rightarrow \Delta_k, \quad \Delta_i \rightarrow -\Delta_j, \quad \Delta_j \rightarrow \Delta_i. \quad (C3) $$

The temporal links live in sites in the corners of the plaquettes and have alternating signs as we go around a plaquette. Upon rotations their positions in the plaquettes will also rotate. This implies that upon rotations the time link variables should change sign. The differences in the temporal direction are unaffected.

$$ A_0 \rightarrow -A_0, \quad \Delta_0 \rightarrow \Delta_0. \quad (C4) $$

By taking all this into consideration it is straightforward to see that the action is invariant under 90° rotations around any axis.

Now we can look into reflection symmetries. While these are easier to visualize, we present a visual representation of reflections on the gauge variables for clarity.

<figure><img src="image_30.png" /><figcaption>FIG. 24: Changes due to a reflection from $k$ to $-k$</figcaption></figure>

The effects of reflection from $k$ to $-k$ on the plaquettes are then

$$ A_{ij} \rightarrow A_{ij}, \quad A_{ik} \rightarrow -A_{ik}, \quad A_{jk} \rightarrow -A_{jk} \quad (C5) $$

and on the differences

$$ \Delta_k \rightarrow -\Delta_k, \quad \Delta_i \rightarrow \Delta_i, \quad \Delta_j \rightarrow \Delta_j. \quad (C6) $$

For this type of transformation both the time links and the temporal differences remain unchanged

$$ A_0 \rightarrow A_0, \quad \Delta_0 \rightarrow \Delta_0. \quad (C7) $$
Again, it is simple to see that the action is invariant for reflections around any direction.

# 1. Correlation Functions

While the action is invariant to reflections and translation, the same is not necessarily true for the correlation functions $\langle L_{ij\tau}(\vec{r}_1)L_{ij\tau}(\vec{r}_2) \rangle$ and $\langle L_{ik\tau}(\vec{r}_1)L_{kj\tau}(\vec{r}_2) \rangle$.

From the earlier results we can deduce that the operators $L_{ij\tau}$ transform in the following way for rotations

$$L_{ij} \xrightarrow{k \to 90^\circ} L_{ij} \qquad L_{ij} \xrightarrow{i \to 90^\circ} L_{ki}^\dagger \qquad L_{ij} \xrightarrow{j \to 90^\circ} L_{jk}. \quad (C8)$$

This then implies that for the correlation operators we have

$$\begin{align}
L_{ij} L_{ij}^\dagger \xrightarrow{k \to 90^\circ} L_{ij} L_{ij}^\dagger \tag{C9} \\
L_{ij} L_{ij}^\dagger \xrightarrow{i \to 90^\circ} L_{ki}^\dagger L_{ki} \xrightarrow{i \to 90^\circ} L_{ij}^\dagger L_{ij} \nonumber \\
L_{ij} L_{ij}^\dagger \xrightarrow{j \to 90^\circ} L_{jk} L_{jk}^\dagger \xrightarrow{j \to 90^\circ} L_{ij}^\dagger L_{ij} \tag{C10}
\end{align}$$

and

$$\begin{align}
L_{ik} L_{kj}^\dagger \xrightarrow{k \to 90^\circ} L_{kj}^\dagger L_{ki}^\dagger \xrightarrow{k \to 90^\circ} L_{ik}^\dagger L_{kj} \tag{C11} \\
L_{ik} L_{kj}^\dagger \xrightarrow{i \to 90^\circ} L_{ij} L_{kj}^\dagger \xrightarrow{i \to 90^\circ} L_{ik}^\dagger L_{kj}^\dagger \xrightarrow{i \to 90^\circ} L_{ij}^\dagger L_{kj}^\dagger \xrightarrow{i \to 90^\circ} L_{ik} L_{kj}^\dagger \tag{C12} \\
L_{ik} L_{kj}^\dagger \xrightarrow{j \to 90^\circ} L_{ik} L_{ij} \xrightarrow{j \to 90^\circ} L_{ik} L_{kj} \xrightarrow{j \to 90^\circ} L_{ik} L_{ij}^\dagger \xrightarrow{j \to 90^\circ} L_{ik} L_{kj}^\dagger \tag{C13}
\end{align}$$

The symmetries present in the theory do not carry on to the correlation functions. For $\langle L_{ij\tau}(\vec{r}_1)L_{ij\tau}(\vec{r}_2) \rangle$ we keep the four-fold symmetry around the $k$ axis, but it breaks down to two-fold symmetry for the $i$ and $j$ axis. The lack of symmetry is even more prevalent for $\langle L_{ik\tau}(\vec{r}_1)L_{kj\tau}(\vec{r}_2) \rangle$, where there is only two-fold symmetry around $k$.

A similar analysis can be done for reflections. The elementary operators $L_{ij\tau}$ transform like

$$L_{ij} \xrightarrow{k \to -k} L_{ij} \qquad L_{ij} \xrightarrow{i \to -i} L_{ij}^\dagger \qquad L_{ij} \xrightarrow{j \to -j} L_{ij}^\dagger. \quad (C14)$$

This then implies that for the correlation operators we have

$$\begin{align}
L_{ij} L_{ij}^\dagger \xrightarrow{k \to -k} L_{ij} L_{ij}^\dagger \qquad L_{ij} L_{ij}^\dagger \xrightarrow{i \to -i} L_{ij}^\dagger L_{ij} \nonumber \\
L_{ij} L_{ij}^\dagger \xrightarrow{j \to -j} L_{ij}^\dagger L_{ij} \tag{C15}
\end{align}$$

and

$$\begin{align}
L_{ik} L_{kj}^\dagger \xrightarrow{k \to -k} L_{ik}^\dagger L_{kj} \qquad L_{ik} L_{kj}^\dagger \xrightarrow{i \to -i} L_{ik}^\dagger L_{kj}^\dagger \nonumber \\
L_{ik} L_{kj}^\dagger \xrightarrow{j \to -j} L_{ik} L_{kj}. \tag{C16}
\end{align}$$

For reflections the symmetries prevail for $\langle L_{ij\tau}(\vec{r}_1)L_{ij\tau}(\vec{r}_2) \rangle$ but, again, the operator $\langle L_{ik\tau}(\vec{r}_1)L_{kj\tau}(\vec{r}_2) \rangle$ proves to have a higher degree of assymetry since the reflection symmetries around $i$ and $j$ break down.

A summary of these findings can be found in the table below.

<table><thead><tr><th>Operator</th><th>Rotation</th><th>Reflection</th></tr></thead><tbody><tr><td>$L_{ij} L_{ij}^\dagger$</td><td>$k^{90^\circ}, i^{180^\circ}$ and $j^{180^\circ}$</td><td>All axis</td></tr><tr><td>$L_{ik} L_{kj}^\dagger$</td><td>$k^{180^\circ}$</td><td>$k \to -k$</td></tr></tbody></table>

[1] K. G. Wilson, Confinement of quarks, Phys. Rev. D 10, 2445 (1974).

[2] A. M. Polyakov, Quark confinement and topology of gauge theories, Nuclear Physics B 120, 429 (1977).

[3] M. Pretko, Generalized electromagnetism of subdimensional particles: A spin liquid story, Physical Review B 96, 10.1103/physrevb.96.035119 (2017).

[4] M. Pretko, Subdimensional particle structure of higher rank u(1) spin liquids, Physical Review B 95, 10.1103/physrevb.95.115139 (2017).

[5] J. Haah, Local stabilizer codes in three dimensions without string logical operators, Phys. Rev. A 83, 042330 (2011).

[6] S. Bravyi and J. Haah, Quantum self-correction in the 3d cubic code model, Phys. Rev. Lett. 111, 200501 (2013).

[7] 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, 10.1103/physrevb.92.035111 (2019).
235136 (2015).

[8] S. Vijay, J. Haah, and L. Fu, Fracton topological order, generalized lattice gauge theory, and duality, Phys. Rev. B 94, 235157 (2016).

[9] R. M. Nandkishore and M. Hermele, Fractons, Annual Review of Condensed Matter Physics 10, 295 (2019).

[10] M. Pretko, X. Chen, and Y. You, Fracton phases of matter, International Journal of Modern Physics A 35, 203003 (2020).

[11] A. Gromov and L. Radzihovsky, Colloquium: Fracton matter, Rev. Mod. Phys. 96, 011001 (2024).

[12] Y. You, Quantum liquids: Emergent higher-rank gauge theory and fractons (2024), arXiv:2403.17074 [cond-mat.str-el].

[13] Z.-C. Gu and X.-G. Wen, A lattice bosonic model as a quantum theory of gravity, arXiv e-prints , gr-qc/0606100 (2006), arXiv:gr-qc/0606100 [gr-qc].

[14] C. Xu, Gapless bosonic excitation without symmetry breaking: An algebraic spin liquid with soft gravitons, Phys. Rev. B 74, 224433 (2006).

[15] C. Xu and C. Wu, Resonating plaquette phases in SU(4) Heisenberg antiferromagnet, Phys. Rev. B 77, 134449 (2008).

[16] A. Rasmussen, Y.-Z. You, and C. Xu, Stable gapless bose liquid phases without any symmetry (2016).

[17] K. Slagle, A. Prem, and M. Pretko, Symmetric tensor gauge theories on curved spaces, Annals of Physics 410, 167910 (2019), arXiv:1807.00827 [cond-mat.str-el].

[18] Y.-H. Du, U. Mehta, D. Nguyen, and D. T. Son, Volume-preserving diffeomorphism as nonabelian higher-rank gauge symmetry, SciPost Physics 12, 050 (2022), arXiv:2103.09826 [cond-mat.str-el].

[19] M. Qi, L. Radzihovsky, and M. Hermele, Fracton phases via exotic higher-form symmetry-breaking, Annals of Physics 424, 168360 (2021), arXiv:2010.02254 [cond-mat.str-el].

[20] N. Seiberg and S.-H. Shao, Exotic U(1) symmetries, duality, and fractons in 3+1-dimensional quantum field theory, SciPost Physics 9, 10.21468/scipostphys.9.4.046 (2020).

[21] P. Gorantla, H. T. Lam, N. Seiberg, and S.-H. Shao, A modified villain formulation of fractons and other exotic theories, Journal of Mathematical Physics 62, 102301 (2021).

[22] P. Gorantla, H. T. Lam, N. Seiberg, and S.-H. Shao, Low-energy limit of some exotic lattice theories and uv/ir mix-

ing, Phys. Rev. B 104, 235116 (2021).

[23] P. Gorantla, H. T. Lam, N. Seiberg, and S.-H. Shao, Global dipole symmetry, compact Lifshitz theory, tensor gauge theory, and fractons, Phys. Rev. B 106, 045112 (2022).

[24] P. Gorantla, H. T. Lam, N. Seiberg, and S.-H. Shao, (2+1)-dimensional compact Lifshitz theory, tensor gauge theory, and fractons, Physical Review B 108, 075106 (2023), arXiv:2209.10030 [cond-mat.str-el].

[25] A. Prem, J. Haah, and R. Nandkishore, Glassy quantum dynamics in translation invariant fracton models, Phys. Rev. B 95, 155133 (2017).

[26] M. Pretko and L. Radzihovsky, Fracton-elasticity duality, Phys. Rev. Lett. 120, 195301 (2018).

[27] C. Xu and C. Wu, Resonating plaquette phases in su(4) heisenberg antiferromagnet, Phys. Rev. B 77, 134449 (2008).

[28] H. Ma, M. Hermele, and X. Chen, Fracton topological order from the higgs and partial-confinement mechanisms of rank-two gauge theory, Physical Review B 98, 10.1103/physrevb.98.035111 (2018).

[29] D. Bulmash and M. Barkeshli, Higgs mechanism in higher-rank symmetric u(1) gauge theories, Physical Review B 97, 10.1103/physrevb.97.235112 (2018).

[30] E. Fradkin and S. H. Shenker, Phase diagrams of lattice gauge theories with higgs fields, Phys. Rev. D 19, 3682 (1979).

[31] D. J. E. Callaway and L. J. Carson, Abelian higgs model: A monte carlo study, Phys. Rev. D 25, 531 (1982).

[32] A. M. Polyakov, Gauge Fields and Strings (Routledge, London, 2017).

[33] A. Prem, S. Vijay, Y.-Z. Chou, M. Pretko, and R. M. Nandkishore, Pinch point singularities of tensor spin liquids, Phys. Rev. B 98, 165140 (2018).

[34] S. Vijay, J. Haah, and L. Fu, Fracton topological order, generalized lattice gauge theory, and duality, Physical Review B 94, 10.1103/physrevb.94.235157 (2016).

[35] A. H. Guth, Existence proof of a nonconfining phase in four-dimensional u(1) lattice gauge theory, Phys. Rev. D 21, 2291 (1980).

[36] K. Osterwalder and E. Seiler, Gauge field theories on a lattice, Annals of Physics 110, 440 (1978).

[37] K. T. K. Chung, R. Flores-Calderón, R. C. Torres, P. Ribeiro, S. Moroz, and P. McClarty, Higgs Phases and Boundary Criticality, arXiv e-prints , arXiv:2404.17001 (2024), arXiv:2404.17001 [cond-mat.str-el].