# Fractonic Berezinskii-Kosterlitz-Thouless transition
from a renormalization group perspective

Kevin T. Grosvenor<EMAIL_ADDRESS>Instituut-Lorentz,
Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands.  Ruben
Lier<EMAIL_ADDRESS>Max Planck Institute for the Physics of Complex
Systems, 01187 Dresden, Germany  Würzburg-Dresden Cluster of Excellence
ct.qmat, Germany  Piotr Surówka<EMAIL_ADDRESS>Department of
Theoretical Physics, Wrocław University of Science and Technology, 50-370
Wrocław, Poland  Institute for Theoretical Physics, University of Amsterdam,
1090 GL Amsterdam, The Netherlands  Dutch Institute for Emergent Phenomena
(DIEP), University of Amsterdam, 1090 GL Amsterdam, The Netherlands

###### Abstract

Proliferation of defects is a mechanism that allows for topological phase
transitions. Such a phase transition is found in two dimensions for the XY-
model, which lies in the Berezinskii-Kosterlitz-Thouless (BKT) universality
class. The transition point can be found using renormalization group analysis.
We apply renormalization group arguments to determine the nature of BKT
transitions for the three-dimensional plaquette-dimer model, which is a model
that exhibits fractonic mobility constraints. We show that an important part
of this analysis demands a modified dimensional analysis that changes the
interpretation of scaling dimensions upon coarse-graining. Using this modified
dimensional analysis we compute the beta functions of the model and predict a
finite critical value above which the fractonic phase melts, proliferating
dipoles. Importantly, the transition point and its value are found
unequivocally within the formalism of renormalization group.

##  I Introduction

Quantum excitations with mobility constraints constitute new phases of matter
called fractons. In the quest of understanding properties of these phases and
their experimental realization it is necessary to develop their macroscopic
description and identify unique properties that distinguish them from other
phases of matter studied in the past. The framework that allows for a
systematic investigation of physical systems at different scales is called the
renormalization group (RG). Unfortunately a direct application of this
framework to fracton models is faced with difficulties due to the so-called
UV/IR mixing phenomenon, which in essence means that, depending on the chosen
direction, low-energy modes can have very high momenta. As a consequence it
was suggested that renormalization group is not applicable to fracton phases
Zhou:2021wsv ; moessneryizhi . However, it was recently argued that this
difficulty can be circumvented by adapting the integration of the high-energy
modes to the symmetries exhibited by the fracton models ethanlake .

Renormalization group analysis of the XY model reveals that it lies in the
universality class of the Berezinskii-Kosterlitz-Thouless (BKT) transition
Kosterlitz_1973 ; berezinskii1971destruction . A distinctive feature of this
universality class is the critical temperature that governs the proliferation
of free topological defects. Applying a similar analysis to fractonic theories
can potentially provide a novel diagnostic of universality classes in theories
with mobility constraints. Our goal is to refine the framework of
renormalization group to study the proliferation of defects in a theory with
mobility constrains. A new universality class of the BKT type was recently
anticipated in the context of superfluids and plaquette-dimer liquids
Zhou:2021wsv .

Dimer models represent lattice systems with degrees of freedom on the links
instead of the nodes dimer1 ; dimer2 ; dimer3 ; dimer4 ; dimer5 ; dimer6 ;
dimer7 ; dimer8 ; dimer9 ; dimer10 . Such models originate from the quest of
understanding magnetic materials and are used to shed light on valence bond
liquids or classical spin ice. In constructing these dimer models one imposes
a constraint on the dimers, namely that each site form a dimer with only one
of its neighbors. Sites that violate this condition are associated with
defects. A site that is not attached to any dimer is called a monomer. A
single monomer cannot move alone, while a pair of monomers between links can
only move along the transverse direction. A generalization of simple dimers
involves trimers and plaquettes. Crucially, a class of plaquette-dimer models
can be mapped to electrostatics with higher-rank tensor electric fields,
considered by Pretko in the context of gapless fractons
pretko_generalized_2017 ; pretko_subdimensional_2017 . In addition the low-
energy effective theory is governed by the physics of defects, i.e., singular
configurations of the fields seiberg_exotic_2020 . The goal of the present
paper is to employ renormalization group analysis to study the transition from
the liquid phase to the ordered configuration. Our formalism can be applied to
a variety of fracton models using the powerful technique of the
renormalization group.

##  II Fractonic plaquette-dimer model

In this work we study the fractonic plaquette-dimer model also studied in Ref.
moessneryizhi . The partition function of this model is given by

$\displaystyle\begin{split}\mathcal{Z}&=\int\mathcal{D}h\,\text{exp}\Bigg{\\{}-\int
d^{3}\mathbf{x}\bigg{[}\frac{\kappa}{2}(\partial_{x}\partial_{y}h)^{2}+\frac{\kappa}{2}(\partial_{z}h)^{2}+2\sum_{I}\alpha_{I}\cos(2\pi
f_{I})\bigg{]}\Bigg{\\}}~{}~{}.\end{split}$ (1)

$f_{I}$ are general functions that enter through a cosine in order to satisfy
the “integer constraint” moessneryizhi  and they have the corresponding
“fugacity” $\alpha_{I}$. In this work we will consider the functions

$\displaystyle
f_{I}=\\{h,a_{x}\partial_{x}h,a_{y}\partial_{y}h\\}~{}~{},~{}~{}I\in\\{0,x,y\\}~{}~{},$
(2)

with $a_{x}$ and $a_{y}$ being the lattice constants. We study the system at
static equilibrium, which is why the time integral does not appear in the
action in Eq. (1), having been canceled by the inverse temperature in
Euclidean signature. Furthermore, the continuum limit was taken in Eq. (1), as
is standard when doing RG for lattice models, such as for example the XY model
herbut_2007 ; tongtong . However, the nature of the model in Eq. (1),
specifically the UV/IR mixing, is such that discreteness will give
quantitative contributions for which there is no limit that can make them
negligible. For simplicity and analytic tractability, the continuum limit will
be used for the qualitative message of this work, but the quantitative effect
of discreteness will also be discussed.

Since the Gaussian term in Eq. (1), proportional to $\kappa$, gives a
vanishing dispersion when $p=0$ or $q=0$, with
$\mathbf{p}=(p,q,k)^{\intercal}$ being the wave-vector, even at low energies
one has to take note of short-wavelength effects. This phenomenon is called
UV/IR mixing uvirmix1 ; uvirmix2 ; you2020fracton ; you2021fractonic , and
because of this effect it is argued in Ref. moessneryizhi  that this model is
beyond the renormalization group paradigm. By considering the correlation
functions, those authors are nevertheless able to derive a critical point
where the low-energy theory is no longer described by a free fractonic theory
but instead by a proliferation of dipoles. Because this proliferation of
defects destroys the quasi-long-range order, it is reminiscent of the BKT
transition Kosterlitz_1973 ; berezinskii1971destruction . Corresponding to
this phase transition is a critical $\kappa_{c}$, above which the coefficients
$\alpha_{x}$ and $\alpha_{y}$ are relevant and grow large in the IR. When this
happens, the field $h$ arranges itself to be in the valley of the cosines with
$I=x,y$, so that at low energies the cosines can be expanded and one is left
with a three-dimensional sine-Gordon model kadanoff  that does not have any
fractonic properties. In this work we show that it is possible to do momentum
shell RG for the model of Eq. (1) by using an RG procedure which is an
extension of the RG procedure in Ref. ethanlake , where RG for the exciton
Bose liquid balentsringexchange ; cenkexu  is considered, which is a quantum
model that also suffers from UV/IR mixing. With this approach the momentum
shells are along the constant energy surface of the fractonic dispersion (see
figure 1). Because of this, one flows towards the axes as opposed to the
origin, and one thus avoids issues related to UV/IR mixing. In sec. III we
discuss this approach and specifically the nature of the dilatation operator
for this RG procedure. We use this in sec. IV to derive the critical point
$\kappa_{c}$. With this approach we are furthermore enabled to compute the
screening effect that the cosine terms have on $\kappa$, i.e., the inverse of
the “fractonic spin-stiffness”, which we find to be absent.

##  III RG procedure

The RG procedure of integrating out high-energy modes in a way which is
adapted to the fractonic plaquette-dimer model was first first proposed and
implemented in ethanlake . However, the crucial subsequent step of rescaling
the low-energy modes back up in order to return to the original form of the
action, but with renormalized parameters, was still done the usual way, which
is _not_ adapted to the fractonic model. In this section, we describe the
appropriate modification to this second step. In essence, what we are pointing
out here is that the concept of _dimensional analysis_ itself is modified near
the fractonic free-field fixed point as compared with the homogeneous one.
This will obviously influence what we mean by _relevant_ , _irrelevant_ , and
_marginal_ in the context of RG.

The most natural way to make manifest the concept in RG of “flowing to the IR”
is to integrate out high-_energy_ shells. These contain only modes with
energies between some high-energy cut-off $\Lambda$ to some slightly lower
energy scale $\Lambda/b$, where $b$ is a number slightly greater than 1. We
emphasize that this is a high-energy shell and not a high-momentum shell, as
it is often called; it simply happens to be the case that high energy and high
momentum are usually equivalent. This is _not_ the case for the fractonic
phase. Nevertheless, the procedure is operationally the same: we integrate out
modes with energy between $\Lambda$ and $\Lambda/b$. In momentum space,
$(p,q,k)$, the constant-energy surface is given by
$\sqrt{(pq)^{2}+k^{2}}=\Lambda$ and depicted in Fig. 1. Whereas the zero-
energy locus is usually just one point, namely the origin in momentum space,
in our case, it is the union of the $p$\- and $q$-axes.

Figure 1: Constant-energy surfaces.

Normally, the constant-energy surfaces are concentric spheres and the RG flow
is towards the origin. Therefore, the RG flow is described by the dilatation
operator in momentum space $D=\mathbf{p}\cdot\nabla_{\mathbf{p}}$. Thus, after
integrating out a shell between energies $\Lambda$ and $\Lambda/b$, we rescale
momenta precisely by the dilatation operator in order to bring $\Lambda/b$
back up to $\Lambda$. In this regard, what is meant by the dimension of an
operator or a parameter is really just its eigenvalue with respect to
dilatation.

In our case, it no longer makes sense to define scaling dimension with respect
to the dilatation operator $D$ since that does _not_ describe the RG flow
towards the IR anymore. Instead, we define a modified dilatation operator
$\tilde{D}$, which describes the flow from one constant-energy surface to
another. First, define the modified momentum vector

$\tilde{\mathbf{p}}=\bigl{(}[p]\,p,[q]\,q,[k]k\bigr{)}^{\intercal}~{}~{},$ (3)

where $[p],[q],[k]$ denote the dimensions of the momentum components. The new
dilatation operator is

$\tilde{D}=\tilde{\mathbf{p}}\cdot\nabla_{\mathbf{p}}~{}~{},$ (4)

under which the eigenvalues of $p$, $q$ and $k$ are, by construction, their
dimensions. We ask that $\tilde{D}$ be orthogonal to the constant-energy
surfaces. Furthermore, we have the freedom to set the dimension of any one
component to 1, which then fixes the remaining two. It is natural to set
$[k]=1$, which fixes the dimensions to be

$\displaystyle[p]$ $\displaystyle=\frac{q^{2}}{p^{2}+q^{2}}~{}~{},$
$\displaystyle[q]$ $\displaystyle=\frac{p^{2}}{p^{2}+q^{2}}~{}~{},$
$\displaystyle[k]$ $\displaystyle=1~{}~{}.$ (5)

This means that if we rescale the energy $\Lambda\rightarrow b\Lambda$, then
the momenta are rescaled according to $p\rightarrow b^{[p]}p$, and similarly
for $q$ and $k$. This is quite different and exotic compared with the standard
dimensional analysis, which would say that $p$ and $q$ both have dimension
$\frac{1}{2}$. Of course, in either case, the dimension of the product $pq$ is
still given by

$\displaystyle[pq]=1.$ (6)

Since the meaning of classically relevant, irrelevant, and marginal depends on
dimensional analysis, the above modified definition of dimensions has some
rather important consequences. One issue is the question of whether gradient
operators without fractonic properties, such as $(\partial_{x}h)^{2}$ and
$(\partial_{y}h)^{2}$, will get generated. The naive form of dimensional
analysis would conclude that these operators are relevant compared with the
operator $(\partial_{x}\partial_{y}h)^{2}$. However, with respect to the
modified dimensional analysis appropriate near the fractonic fixed point,
$(\partial_{x}h)^{2}$ and $(\partial_{y}h)^{2}$ do not even have constant
dimensions and, in fact, their highest dimension is 2, which is the same as
the dimension of the operator $(\partial_{x}\partial_{y}h)^{2}$. This may seem
to suggest that these operators are no more relevant than the fractonic one.
However, consider, for example, the dimension of the operator
$(\partial_{y}h)^{2}$ in the region near the $q$-axis, which is where it could
destroy UV/IR mixing. Here, this dimension vanishes and from the point of view
of the potential destruction of UV/IR mixing, this operator is even more
relevant than in the naive dimensional analysis. We conclude that the
fractonic theory is highly susceptible to being destroyed by simple gradient
terms without fractonic properties. We thus have to start with a pristine
fractonic theory without the ordinary gradient terms. In this case,
dimensional analysis alone cannot determine whether or not the ordinary
gradient terms do in fact get generated. We then ask what kind of theory the
cosine perturbations induce in the IR. To answer this question, we will now
perform the RG analysis along the lines we described above.

##  IV Fractonic plaquette-dimer melting

With the RG procedure described in sec. III, we can now derive the critical
point where the fractonic dimer-plaquette liquid becomes unstable, which was
discussed in Ref. moessneryizhi . This critical point is due to the cosine
operators in Eq. (1) with $I=x$ and $I=y$ becoming relevant. In the
appendices, following Refs. Herbut_2003 ; seradjeh ; herbut_2007 , we consider
the effect of a general cosine term on the renormalization of the coefficients
when one splits the field $h$ into

$\displaystyle h=h_{L}+h_{H}~{}~{},$ (7)

where $h_{H}$ corresponds to the high-energy modes in the momentum shell that
are being integrated out, leading to a renormalization of the coefficients for
the low-energy theory with modes $h_{L}$. Taking $I=x$, one finds that the
fugacity coefficient $\alpha_{x}$ experiences the following renormalization:

$\displaystyle\alpha_{x}(b)=b^{2}\alpha_{x}e^{-\frac{1}{2}g_{H}^{(xx)}(0)}~{}~{},$
(8)

where $g^{(IJ)}_{H}(0)$ is the correlator of the high-energy modes for
operators $f_{I}$ and $f_{J}$. The factor $b^{2}$ is consistent with the
rescaling for the momenta described in sec. III, specifically eqs. (5) and
(6). We thus must compute

$\displaystyle g^{(xx)}_{H}(0)$
$\displaystyle=\frac{a_{x}^{2}}{\kappa}\int_{\Lambda\rightarrow\Lambda/b}\frac{d^{3}\mathbf{p}}{2\pi}\frac{p^{2}}{(pq)^{2}+k^{2}}~{}~{}.$
(9)

To do this, we define $\tilde{q}=q|p|$ and $\tilde{k}=k$ and introduce polar
coordinates with the following radial coordinate:

$\displaystyle\tilde{q}^{\prime 2}\equiv\tilde{q}^{2}+\tilde{k}^{2}~{}~{}.$
(10)

With this coordinate transformation the momentum shell corresponding to the RG
procedure of this work goes from $\tilde{q}^{\prime}=\Lambda$ to
$\tilde{q}^{\prime}=\Lambda/b$. To make sure that the momentum shells of the
cover an as wide as possible angular range we impose that
$\tilde{\Lambda}=a_{x}a_{y}\Lambda/\pi^{2}\ll 1$. This is analogous to
requiring that a spherical momentum shell lies completely inside of the box
formed by the momentum cutoff. Integrating over $\tilde{q}^{\prime}$ and $p$,
we then find111In Appendix C, the integral in Eq. (9) is computed without
assuming $\tilde{\Lambda}\ll 1$. The final result Eq. (11) is retrieved in the
limit $\tilde{\Lambda}\rightarrow 0$.

$\displaystyle\begin{split}g^{(xx)}_{H}(0)&=\frac{a_{x}^{2}}{\kappa}\int^{\pi/a_{x}}_{-\pi/a_{x}}dp\int^{\Lambda}_{\Lambda/b}d\tilde{q}^{\prime}\frac{p^{2}}{|p|\tilde{q}^{\prime}}=\frac{\pi^{2}}{\kappa}\log(b)~{}~{}.\end{split}$
(11)

The fact that the integral in Eq. (11) depends on the UV cutoff $a_{x}$ even
though this RG flow may well happen deeper into the IR is characteristic of
fractonic models where there is UV/IR mixing ethanlake . For $\alpha_{x}$ we
thus get the RG flow:

$\displaystyle\frac{\partial\alpha_{x}}{\partial\log(b)}=\biggl{(}2-\frac{\pi^{2}}{2\kappa}\biggr{)}\alpha_{x}~{}~{}.$
(12)

The same flow occurs for $\alpha_{y}$. From Eq. (12) we learn that the
critical point for this BKT-like transition is given by

$\displaystyle\kappa_{c}=\pi^{2}/4~{}~{}.$ (13)

When formulating the model in Eq. (1), the continuum assumption was used. This
limit works only when the lattice constant can be assumed to always be
unimportant for the RG computations. However, due to the UV/IR mixing, the IR
theory contains momentum modes near the axes that reach the momentum cutoff
set by the inverse of the lattice constant, so that this assumption is never
valid. Because of this, incorporating lattice effects will modify Eq. (13) so
that a different $\kappa^{(L)}_{c}$ is found. Specifically, the result found
in Appendix D is given by

$\displaystyle\kappa^{(L)}_{c}=2~{}~{},$ (14)

so that the lattice effects reduce the critical parameter by approximately
$20\%$ relative to the continuous limit.

We now consider the flow for $\alpha_{0}$:

$\displaystyle\alpha_{0}(b)=b^{2}\alpha_{0}e^{-\frac{1}{2}g_{H}^{(00)}(0)}~{}~{}.$
(15)

$g_{H}^{(00)}(0)$ is more difficult to compute than $g_{H}^{(xx)}(0)$ because
of the absence of the $p^{2}$ term in the denominator, which ensured that the
integrand in Eq. (9) vanished for $p\rightarrow 0$. In Appendix B, we work out
Eq. (15) to find the flow

$\displaystyle\frac{\partial\alpha_{0}}{\partial\log(b)}=\left[2-\frac{1}{\kappa}\log\left(\frac{2}{\tilde{\Lambda}}\right)\right]\alpha_{0}~{}~{}.$
(16)

Because $\tilde{\Lambda}\ll 1$ we find that the coefficient $\alpha_{0}$ can
always be considered to be irrelevant, confirming the result in Ref.
moessneryizhi . To summarize, we have derived a similar prediction as that of
Ref. moessneryizhi , namely that there is a finite critical value $\kappa_{c}$
above which the fractonic phase melts via dipole proliferation.We find only a
difference in the precise value of $\kappa_{c}$. Crucially, we have
established this result now firmly within the formalism of the renormalization
group by working with momentum shells that adhere to UV/IR mixing.

However, now that we know that the $I=0$ cosine term is irrelevant, one may
still worry that while this irrelevant term flows to the IR shrinking rapidly,
it will act as a dangerously irrelevant operator by giving rise to a new term
in the Gaussian part in Eq. (1). Let us imagine that this Gaussian part
$S_{G}$ gets modified so that it now includes the following $\gamma$-term:

$\displaystyle\begin{split}S^{\prime}_{G}&=\frac{\kappa}{2}\left[(\partial_{x}\partial_{y}h)^{2}+(\partial_{z}h)^{2}\right]+\frac{\gamma}{2}\left[(\partial_{x}h)^{2}+(\partial_{y}h)^{2}\right]~{}~{},\end{split}$
(17)

Even for small $\gamma$, this term violates the fractonic properties that were
there for $\gamma=0$ and thus undoes the arguments related to the stability of
the fractonic plaquette-dimer liquid as these arguments relied on UV/IR
mixing. We can extract for $I=0$, i.e. $f_{I}=h$, the integral that gives the
generation of an initially zero $\gamma$ due to integrating out a momentum
shell. This integral is given by

$\displaystyle\gamma(b)$ $\displaystyle=(2\pi)^{2}\alpha_{0}^{2}(b)\int
d^{3}\mathbf{x}\,(x^{2}+y^{2})(\exp(g^{(hh)}_{H}(\mathbf{x}))-1)$
$\displaystyle\approx 2\pi\frac{\alpha_{0}^{2}(b)}{\kappa}\int
d^{3}\mathbf{x}\int_{\Lambda\rightarrow\Lambda/b}d^{3}\mathbf{p}\frac{(x^{2}+y^{2})\exp(i\mathbf{p}\cdot\mathbf{x})}{(pq)^{2}+k^{2}}~{}~{}.$
(18)

This integral is complicated due to the nature of the cut-off surfaces.
Therefore, we deform the region of integration in momentum space such that $k$
is integrated out first from $-\infty$ to $+\infty$. We could have employed
this same trick on the previous integrals as well; this does not affect the
beta function. After integrating both $k$ and $z$, we find

$\displaystyle\begin{split}&\gamma(b)=\frac{(2\pi)^{2}\alpha_{0}^{2}(b)}{\kappa}\int
dx\,dy\int_{\Lambda\rightarrow\Lambda/b}dp\,dq\,\frac{(x^{2}+y^{2})e^{ipx+iqy}}{(pq)^{2}}~{}~{}.\end{split}$
(19)

We move to polar coordinates both for real space as well as for momentum space
and first integrate over the angle $\theta$ in real space:

$\displaystyle\begin{split}&\gamma(b)=\frac{(2\pi)^{3}\alpha_{0}^{2}(b)}{\kappa}\int_{0}^{\infty}dr\int_{\Lambda\rightarrow\Lambda/b}d\rho
d\phi\frac{r^{3}J_{0}(r\rho)}{\frac{1}{4}\sin^{2}(2\phi)\rho^{3}}~{}~{}.\end{split}$
(20)

Note that $\sin^{2}(2\phi)>0$, because the angular integral is cut off where
the shell hits the box formed by the momentum cutoff. We then change the
integration variables as $\tilde{\rho}=\frac{1}{2}|\sin(2\phi)|\rho$ and
$\tilde{r}=2r/|\sin(2\phi)|$. The integral over $\phi$ is finite, being
regulated by the inverse lattice constants and the cut-off surface, and
positive. Thus, we simply denote the $\phi$ integral as
$I(a_{x},a_{y},\Lambda)$:

$\displaystyle\begin{split}&\gamma(b)=\frac{(2\pi)^{3}\alpha_{0}^{2}(b)I(a_{x},a_{y},\Lambda)}{\kappa}\int_{0}^{\infty}d\tilde{r}\int^{\Lambda}_{\Lambda/b}d\tilde{\rho}\frac{\tilde{r}^{3}J_{0}(\tilde{r}\tilde{\rho})}{\tilde{\rho}^{3}}=0~{}~{}.\end{split}$
(21)

The integral in Eq. (21) is divergent, but can be regularized by smoothing the
momentum shell as follows:

$\displaystyle\int^{\Lambda}_{\Lambda/b}d\tilde{\rho}\rightarrow\int^{\infty}_{0}d\tilde{\rho}\left(\frac{\Lambda^{2n}}{\tilde{\rho}^{2n}+\Lambda^{2n}}-\frac{(\Lambda/b)^{2n}}{\tilde{\rho}^{2n}+(\Lambda/b)^{2n}}\right)~{}~{},$
(22)

where for theories with quadratic correlators one takes $n=1$ Herbut_2003 ;
herbut_2007 ; seradjeh . The lowest value for $n$ that regularizes Eq. (21) is
$n=2$. Taking $n=2$ yields

$\displaystyle\gamma(b)=0+\mathcal{O}\left((\log(b))^{2}\right)~{}~{}.$ (23)

We thus conclude that the fractonic UV properties of the model are not
destroyed by a dangerously irrelevant cosine operator when flowing to the IR.
Note that the vanishing integral that leads to Eq. (23) is not specific to
fractonic theories, because the vanishing relies only on the quartic nature of
the correlation function that remains after the third dimension is integrated
out in Eq. (19). This means that for two-dimensional theories with quartic
dispersion, such as the vector sine-Gordon model Zhai_2019 , which describes
dislocation-induced melting of two-dimensional crystals Kosterlitz_1973 ;
kthny1 ; kthny2 ; kthny3 , a similar argument can preclude dangerously
irrelevant cosine operators.

##  V Conclusion

In this work, we looked at a fractonic BKT-transition in three dimensions that
was first mentioned in Ref. moessneryizhi , where a fractonic plaquette-dimer
liquid with algebraic correlations melts into a disordered phase without
fractonic properties. With an RG scheme that is an extension of the one
considered in Ref. ethanlake , we showed that one can derive the critical
properties of this model. This is done through momentum shell RG by
considering the renormalization of the cosine perturbations coming from the
contributions linear in the fugacity of the corresponding cosine operator. We
also considered the effect of contributions second order in the fugacity.
Specifically, we considered the possibility that the simple cosine operator
acts as a dangerously irrelevant operator that destroys the fractonic
dispersion, thereby removing the UV/IR mixing nature of the model in the IR.
Within the scope of our simplified computation, we found that this does not
happen and for this reason we do not find anything preventing the fractonic
properties to live on in the IR. The fractonic properties are then only
destroyed when $\kappa>\kappa_{c}$, in which case a fractonic BKT transition
occurs that causes a flow towards a disordered phase.

_Acknowledgments.–_ We acknowledge discussions with Renato Dantas, Roderich
Moessner, Andrés Schlief Raether, Vijay Shenoy, and Yizhi You. We especially
thank Ethan Lake for explaining his work to us. KTG has received funding from
the European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement No 101024967. RL was supported, in
part, by the cluster of excellence ct.qmat (EXC 2147, project-id 390858490).
PS acknowledges the support of the NCN Sonata Bis grant 2019/34/E/ST3/00405
and the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) Klein
grant via NWA route 2.

##  Appendix A Low-energy effective action

In our RG computations, we will closely follow Ref. herbut_2007  (see also
Refs. Herbut_2003 ; seradjeh ), which performs the same sort of computations
but for the standard BKT transition.

We start with the partition function of our model expanded near the fractonic
free-field fixed point:

$\displaystyle\mathcal{Z}$ $\displaystyle=\int\mathcal{D}h\,e^{-S},$
$\displaystyle S$ $\displaystyle=K+V,$ (24)

where $S$ is the static action and $K$ and $V$ are the kinetic and potential
energies,

$\displaystyle K$ $\displaystyle=\frac{\kappa}{2}\int
d^{3}\mathbf{x}\,\bigl{[}\bigl{(}\partial_{x}\partial_{y}h\bigr{)}^{2}+\bigl{(}\partial_{z}h\bigr{)}^{2}\bigr{]},$
$\displaystyle V$ $\displaystyle=2\int
d^{3}\mathbf{x}\sum_{I}\alpha_{I}\cos(2\pi f_{I}),$ (25)

where the label $I$ takes on values $I=\\{0,x,y\\}$ and where

$\displaystyle f_{0}$ $\displaystyle=h~{}~{},$ $\displaystyle f_{x}$
$\displaystyle=a_{x}\partial_{x}h~{}~{},$ $\displaystyle f_{y}$
$\displaystyle=a_{y}\partial_{y}h~{}~{}.$ (26)

First, let us expand out the contribution of potential to the partition
function:

$\displaystyle\mathcal{Z}$
$\displaystyle=\int\mathcal{D}h\,e^{-K}\sum_{n=0}^{\infty}\frac{V^{n}}{n!}~{}~{},$
$\displaystyle V^{n}$ $\displaystyle=2^{n}\prod_{j=0}^{n}\int
d^{3}\mathbf{x}_{j}\sum_{I_{j}}\alpha_{I_{j}}\cos\bigl(2\pi
f_{I_{j}}\bigr{missing})~{}~{}.$ (27)

We will only keep up to order $n=2$, so we can write these terms out
explicitly:

$\displaystyle V^{2}$ $\displaystyle=4\sum_{I,J}\alpha_{I}\alpha_{J}\int
d^{3}\mathbf{x}\,d^{3}\mathbf{x}^{\prime}\,\cos[2\pi
f_{I}(x)\bigr{]}\,\cos[2\pi f_{J}(x^{\prime})\bigr{]}~{}~{}.$ (28)

To integrate out the momentum shells, as usual, we divide the field $h$ into
high- and low-energy parts:

$\displaystyle h=h_{L}+h_{H}~{}~{}.$ (29)

The modes $h_{H}$ have energies within some energy shell between energies
$\Lambda/b$ and $\Lambda$, where $b$ is a number slightly greater than 1. We
integrate these modes out leaving behind just the low-energy modes $h_{L}$.
The remaining momenta have to then be rescaled back up from $\Lambda/b$ to
$\Lambda$. The latter step similarly scales the couplings according to their
scaling dimension. As we have argued in the main text, the scaling that is
appropriate near the fractonic fixed point is a type of dilatation but towards
the $x$ and $y$ axes, rather than towards the origin.

Since $h_{L}$ and $h_{H}$ have no overlap there are no cross terms between
them arising from the free kinetic part of the action. However, we must still
expand the cosine terms:

$\displaystyle\cos(2\pi f_{I})$
$\displaystyle=\cos[2\pi\bigl{(}f_{I}^{L}+f_{I}^{H}\bigr{)}\bigr{]}=\cos(2\pi
f_{I}^{L})\,\cos(2\pi f_{I}^{H})-\sin(2\pi f_{I}^{L})\,\sin(2\pi
f_{I}^{H})~{}~{}.$ (30)

The partition function now reads

$\displaystyle\mathcal{Z}$
$\displaystyle=\int\mathcal{D}h_{L}\,e^{-K_{L}}\int\mathcal{D}h_{H}\,e^{-K_{H}}\biggl{\\{}1+2\sum_{I}\alpha_{I}\int
d^{3}\mathbf{x}\,\bigl{[}\cos(2\pi f_{I}^{L})\,\cos(2\pi f_{I}^{H})-\sin(2\pi
f_{I}^{L})\,\sin(2\pi f_{I}^{H})\bigr{]}$
$\displaystyle\quad+2\sum_{I,J}\alpha_{I}\alpha_{J}\int
d^{3}\mathbf{x}\,d^{3}\mathbf{x}^{\prime}\,\bigl{[}\cos(2\pi
f_{I}^{L})\,\cos(2\pi f_{I}^{H})-\sin(2\pi f_{I}^{L})\,\sin(2\pi
f_{I}^{H})\bigr{]}(x)$ $\displaystyle\hskip 106.69783pt\times\bigl{[}\cos(2\pi
f_{J}^{L})\,\cos(2\pi f_{J}^{H})-\sin(2\pi f_{J}^{L})\,\sin(2\pi
f_{J}^{H})\bigr{]}(x^{\prime})+O(V^{3})\biggr{\\}}~{}~{}.$ (31)

Integrating out $h_{H}$ amounts to taking the expectation value with respect
to $K_{H}$, which we denote by

$\langle\mathcal{O}\rangle_{H}=\frac{\int\mathcal{D}h_{H}\,e^{-K_{H}}\,\mathcal{O}}{\int\mathcal{D}h_{H}\,e^{-K_{H}}}~{}~{},$
(32)

for an arbitrary functional $\mathcal{O}$ of $h_{H}$.

Since there is no tadpole for $h_{H}$ (no term linear in $h_{H}$), the
expectation value of any function which is odd in $h_{H}$ vanishes. In other
words,

$\displaystyle\bigl{\langle}\sin[2\pi
f_{I}^{H}(x)\bigr{]}\bigr{\rangle}_{H}=\bigl{\langle}\sin[2\pi
f_{I}^{H}(x)\bigr{]}\cos[2\pi
f_{J}^{H}(x^{\prime})\bigr{]}\bigr{\rangle}_{H}=0~{}~{}.$ (33)

Therefore, the only terms which survive in the partition function are

$\displaystyle\mathcal{Z}_{L}$
$\displaystyle=\int\mathcal{D}h_{L}\,e^{-K_{L}}\biggl{\\{}1+2\sum_{I}\alpha_{I}\int
d^{3}\mathbf{x}\,\cos\bigl(2\pi
f_{I}^{L}\bigr{missing})\bigl{\langle}\cos\bigl(2\pi
f_{I}^{H}\bigr{missing})\bigr{\rangle}_{H}$
$\displaystyle\quad+2\sum_{I,J}\alpha_{I}\alpha_{J}\int
d^{3}\mathbf{x}\,d^{3}\mathbf{x}^{\prime}\,\cos[2\pi
f_{I}^{L}(x)\bigr{]}\cos[2\pi
f_{J}^{L}(x^{\prime})\bigr{]}\bigl{\langle}\cos[2\pi
f_{I}^{H}(x)\bigr{]}\cos[2\pi f_{J}^{H}(x^{\prime})\bigr{]}\bigr{\rangle}_{H}$
$\displaystyle\quad+2\sum_{I,J}\alpha_{I}\alpha_{J}\int
d^{3}\mathbf{x}\,d^{3}\mathbf{x}^{\prime}\,\sin[2\pi
f_{I}^{L}(x)\bigr{]}\sin[2\pi
f_{J}^{L}(x^{\prime})\bigr{]}\bigl{\langle}\sin[2\pi
f_{I}^{H}(x)\bigr{]}\sin[2\pi
f_{J}^{H}(x^{\prime})\bigr{]}\bigr{\rangle}_{H}+O(V^{3})\biggr{\\}}~{}~{}.$
(34)

The expectation values of these trigonometric functions are related to the
basic 2-point function

$g_{H}^{(IJ)}(x)=(2\pi)^{2}\bigl{\langle}f_{I}^{H}(x)\,f_{J}^{H}(0)\bigr{\rangle}_{H}~{}~{},$
(35)

via Wick’s theorem:

$\displaystyle\bigl{\langle}\cos[2\pi f_{I}^{H}(x)\bigr{]}\bigr{\rangle}_{H}$
$\displaystyle=e^{-\frac{1}{2}g_{H}^{(II)}(0)}~{}~{},$ (36a)
$\displaystyle\bigl{\langle}\cos[2\pi f_{I}^{H}(x)\bigr{]}\cos[2\pi
f_{J}^{H}(x^{\prime})\bigr{]}\bigr{\rangle}_{H}$
$\displaystyle=e^{-\frac{1}{2}\bigl{[}g_{H}^{(II)}(0)+g_{H}^{(JJ)}(0)\bigr{]}}\cosh[g_{H}^{(IJ)}(x-x^{\prime})\bigr{]}~{}~{},$
(36b) $\displaystyle\bigl{\langle}\sin[2\pi f_{I}^{H}(x)\bigr{]}\sin[2\pi
f_{J}^{H}(x^{\prime})\bigr{]}\bigr{\rangle}_{H}$
$\displaystyle=e^{-\frac{1}{2}\bigl{[}g_{H}^{(II)}(0)+g_{H}^{(JJ)}(0)\bigr{]}}\sinh\bigl{[}g_{H}^{(IJ)}(x-x^{\prime})\bigr{]}~{}~{}.$
(36c)

Therefore, the low-energy partition function reads

$\displaystyle\mathcal{Z}_{L}$
$\displaystyle=\int\mathcal{D}h_{L}\,e^{-K_{L}}\biggl{\\{}1+2\sum_{I}\alpha_{I}e^{-\frac{1}{2}g_{H}^{(II)}(0)}\int
d^{3}\mathbf{x}\,\cos\bigl(2\pi f_{I}^{L}\bigr{missing})$
$\displaystyle\quad+2\sum_{I,J}\alpha_{I}\alpha_{J}e^{-g_{H}^{(IJ)}(0)}\int
d^{3}\mathbf{x}\,d^{3}\mathbf{x}^{\prime}\,\cos[2\pi
f_{I}^{L}(x)\bigr{]}\cos[2\pi
f_{J}^{L}(x^{\prime})\bigr{]}\cosh[g_{H}^{(IJ)}(x-x^{\prime})\bigr{]}$
$\displaystyle\quad+2\sum_{I,J}\alpha_{I}\alpha_{J}e^{-g_{H}^{(IJ)}(0)}\int
d^{3}\mathbf{x}\,d^{3}\mathbf{x}^{\prime}\,\sin[2\pi
f_{I}^{L}(x)\bigr{]}\sin[2\pi
f_{J}^{L}(x^{\prime})\bigr{]}\sinh\bigl{[}g_{H}^{(IJ)}(x-x^{\prime})\bigr{]}+O(V^{3})\biggr{\\}}~{}~{}.$
(37)

We can re-exponentiate this to define a low-energy potential term,

$\displaystyle V_{L}$ $\displaystyle=2\int
d^{3}\mathbf{x}\sum_{I}\alpha_{I}e^{-\frac{1}{2}g_{H}^{(II)}(0)}\cos\bigl(2\pi
f_{I}^{L}\bigr{missing})$ $\displaystyle\quad+2\int
d^{3}\mathbf{x}\,d^{3}\mathbf{x}^{\prime}\sum_{I,J}\alpha_{I}\alpha_{J}e^{-\frac{1}{2}\bigl{[}g_{H}^{(II)}(0)+g_{H}^{(JJ)}(0)\bigr{]}}\biggl{\\{}\cos[2\pi
f_{I}^{L}(x)\bigr{]}\cos[2\pi
f_{J}^{L}(x^{\prime})\bigr{]}\cosh[g_{H}^{(IJ)}(x-x^{\prime})\bigr{]}$
$\displaystyle\hskip 28.45274pt+\sin[2\pi f_{I}^{L}(x)\bigr{]}\sin[2\pi
f_{J}^{L}(x^{\prime})\bigr{]}\sinh\bigl{[}g_{H}^{(IJ)}(x-x^{\prime})\bigr{]}-\cos[2\pi
f_{I}^{L}(x)\bigr{]}\,\cos[2\pi
f_{J}^{L}(x^{\prime})\bigr{]}\biggr{\\}}+\cdots~{}~{},$ (38)

where $\cdots$ are higher-order terms involving more factors of $\alpha_{I}$.
We can massage this into the form

$\displaystyle V_{L}$ $\displaystyle=2\int
d^{3}\mathbf{x}\sum_{I}\alpha_{I}e^{-\frac{1}{2}g_{H}^{(II)}(0)}\cos\bigl(2\pi
f_{I}^{L}\bigr{missing})$ $\displaystyle\quad+2\int
d^{3}\mathbf{x}\,d^{3}\mathbf{x}^{\prime}\sum_{I,J}\alpha_{I}\alpha_{J}e^{-\frac{1}{2}\bigl{[}g_{H}^{(II)}(0)+g_{H}^{(JJ)}(0)\bigr{]}}\biggl{\\{}\cos[2\pi\bigl{(}f_{I}^{L}(x)+f_{J}^{L}(x^{\prime})\bigr{)}\bigr{]}\Bigl{(}e^{-g_{H}^{(IJ)}(x-x^{\prime})}-1\Bigr{)}$
$\displaystyle\hskip
199.16928pt+\cos[2\pi\bigl{(}f_{I}^{L}(x)-f_{J}^{L}(x^{\prime})\bigr{)}\bigr{]}\Bigl{(}e^{g_{H}^{(IJ)}(x-x^{\prime})}-1\Bigr{)}\biggr{\\}}~{}~{}.$
(39)

Finally, we rescale momenta by a factor of $b$ in order to bring the cut-off
back up from $\Lambda/b$ to $\Lambda$ and be able to compare the renormalized
action to the original one. Each factor of $\alpha_{I}$ gets multiplied by
$b^{2}$ since they have dimension 2, as argued in the main text. Let $V(b)$ be
the renormalized potential as a function of $b$:

$\displaystyle V(b)$ $\displaystyle=2\int
d^{3}\mathbf{x}\sum_{I}\alpha_{I}\,b^{2}e^{-\frac{1}{2}g_{H}^{(II)}(0)}\cos\bigl(2\pi
f_{I}^{L}\bigr{missing})$ $\displaystyle\quad+2\int
d^{3}\mathbf{x}\,d^{3}\mathbf{x}^{\prime}\sum_{I,J}\alpha_{I}\alpha_{J}\,b^{4}e^{-\frac{1}{2}\bigl{[}g_{H}^{(II)}(0)+g_{H}^{(JJ)}(0)\bigr{]}}\biggl{\\{}\cos[2\pi\bigl{(}f_{I}^{L}(x)+f_{J}^{L}(x^{\prime})\bigr{)}\bigr{]}\Bigl{(}e^{-g_{H}^{(IJ)}(x-x^{\prime})}-1\Bigr{)}$
$\displaystyle\hskip
210.55022pt+\cos[2\pi\bigl{(}f_{I}^{L}(x)-f_{J}^{L}(x^{\prime})\bigr{)}\bigr{]}\Bigl{(}e^{g_{H}^{(IJ)}(x-x^{\prime})}-1\Bigr{)}~{}~{}.$
(40)

The renormalization of the original cosine terms in the potential are
retrieved from the first line above and we will focus on these terms in the
following sections. The renormalized fugacities read

$\alpha_{I}(b)=\alpha_{I}b^{2}e^{-\frac{1}{2}g_{H}^{(II)}(0)}~{}~{}.$ (41)

This is an implicit equation that requires the evaluation of $\alpha_{I}$ at
some scale. As usual, what we are interested in is the beta function, which is
the logarithmic derivative with respect to $b$.222Our definition of the beta
function is negative of the standard one because by some historical quirk the
usual definition of the beta function flows towards the UV.

Importantly, however, additional terms are generated in the subsequent lines
of $V(b)$ and we will clearly have to address these as well. We will do that
afterwards.

##  Appendix B Renormalization of the simple cosine term

For the simple cosine term, which is when $I=0$ and $f_{0}=h$, we have

$\displaystyle g_{H}^{(00)}(0)=\lim_{x^{\prime}\rightarrow
x}(2\pi)^{2}\bigl{\langle}h(x)h(x^{\prime})\bigr{\rangle}=(2\pi)^{2}\int\frac{d^{3}\mathbf{p}}{(2\pi)^{3}}\frac{1}{\kappa\bigl{[}(pq)^{2}+k^{2}\bigr{]}}~{}~{},$
(42)

where $\mathbf{p}=(p,q,k)$ is the momentum vector in three dimensions. Let us
pass to dimensionless variables by defining

$\displaystyle\tilde{p}$ $\displaystyle=\frac{a_{x}}{\pi}\,p~{}~{},$
$\displaystyle\tilde{q}$ $\displaystyle=\frac{a_{y}}{\pi}\,q~{}~{},$
$\displaystyle\tilde{k}$ $\displaystyle=\frac{a_{x}a_{y}}{\pi^{2}}\,k~{}~{}.$
(43)

In this way, the Brillouin zone in the $pq$-plane is translated into the
square $|\tilde{p}|,|\tilde{q}|\leq 1$. The expression of the integral
actually remains the same, just with tildes on the variables:

$\displaystyle
g_{H}^{(00)}(0)=\frac{1}{2\pi\kappa}\int\frac{d^{3}\tilde{\mathbf{p}}}{(\tilde{p}\tilde{q})^{2}+\tilde{k}^{2}}~{}~{}.$
(44)

The shell goes from $\sqrt{(pq)^{2}+k^{2}}=\frac{\Lambda}{b}$ to $\Lambda$,
which, in these dimensionless variables corresponds to333In fact, it is much
easier to simply integrate $k$ from $-\infty$ to $+\infty$ to begin with. This
effectively reduces the problem to two dimensions. Though this is technically
not in keeping with integrating out modes between constant-energy surfaces, it
turns out to give the same result for the beta function of $\alpha_{I}$.

$\displaystyle\frac{\tilde{\Lambda}}{b}\leq\sqrt{(\tilde{p}\tilde{q})^{2}+\tilde{k}^{2}}\leq\tilde{\Lambda}~{}~{},$
$\displaystyle\tilde{\Lambda}$ $\displaystyle=\frac{\Lambda
a_{x}a_{y}}{\pi^{2}}~{}~{}.$ (45)

But, we must also restrict to the Brillouin zone in the $pq$-plane. Change
variables from $(\tilde{p},\tilde{q},\tilde{k})$ to
$(\tilde{p},\tilde{k},\tilde{\ell})$, where $\tilde{\ell}$ is defined as

$\tilde{\ell}=\tilde{p}\tilde{q}~{}~{}.$ (46)

Remembering to include the determinant of the Jacobian for this change of
variables, which is $\frac{1}{|\tilde{p}|}$, the integral becomes

$\displaystyle
g_{H}^{(00)}(0)=\frac{1}{2\pi\kappa}\int\frac{d\tilde{p}}{|\tilde{p}|}\frac{d\tilde{k}\,d\tilde{\ell}}{\tilde{k}^{2}+\tilde{\ell}^{2}}=\frac{4}{\pi\kappa}\int_{>0}\frac{d\tilde{p}}{\tilde{p}}\frac{d\tilde{k}\,d\tilde{\ell}}{\tilde{k}^{2}+\tilde{\ell}^{2}}~{}~{},$
(47)

where the subscript $>0$ on the integral means restrict to the positive
octant.

The limits of integration in the shell for the $\tilde{p}$ variable is
$\tilde{\ell}$ to 1. Thus,

$\displaystyle
g_{H}^{(00)}(0)=\frac{4}{\pi\kappa}\int_{>0}\frac{d\tilde{k}\,d\tilde{\ell}}{\tilde{k}^{2}+\tilde{\ell}^{2}}\int_{\tilde{\ell}}^{1}\frac{d\tilde{p}}{\tilde{p}}=\frac{4}{\pi\kappa}\int_{>0}\frac{d\tilde{k}\,d\tilde{\ell}}{\tilde{k}^{2}+\tilde{\ell}^{2}}\,\log\frac{1}{\tilde{\ell}}~{}~{}.$
(48)

Switch to polar coordinates

$\displaystyle\tilde{k}$ $\displaystyle=\tilde{\lambda}\,\cos\theta~{}~{},$
$\displaystyle\tilde{\ell}$ $\displaystyle=\tilde{\lambda}\,\sin\theta~{}~{},$
(49)

where, since we are in the positive quadrant, $0\leq\theta\leq\frac{\pi}{2}$.
Then,

$\displaystyle g_{H}^{(00)}(0)$
$\displaystyle=\frac{4}{\pi\kappa}\int_{\tilde{\Lambda}/b}^{\tilde{\Lambda}}\frac{d\tilde{\lambda}}{\tilde{\lambda}}\int_{0}^{\pi/2}d\theta\,\log\biggl(\frac{1}{\tilde{\lambda}\,\sin\theta}\biggr{missing})$
$\displaystyle=\frac{2}{\kappa}\int_{\tilde{\Lambda}/b}^{\tilde{\Lambda}}\frac{d\tilde{\lambda}}{\tilde{\lambda}}\,\log\biggl(\frac{2}{\tilde{\lambda}}\biggr{missing})$
$\displaystyle=\frac{2}{\kappa}\biggl{[}\log\biggl(\frac{2}{\tilde{\Lambda}}\biggr{missing})+\frac{1}{2}\log(b)\biggr{]}\log(b)~{}~{}.$
(50)

Therefore, the renormalized fugacity is

$\displaystyle\alpha_{0}(b)$
$\displaystyle=\alpha_{0}\,b^{2-\frac{1}{\kappa}\bigl{[}\log\bigl(\frac{2}{\tilde{\Lambda}}\bigr{missing})+\frac{1}{2}\log(b)\bigr{]}}~{}~{},$
(51)

and the logarithmic derivative at $b=1$ evaluates to444Unlike in relativistic
theories in which multilogs appear only at two and higher loops, here a
multilog has appeared already at one loop. We treat this in precisely the same
way as in nonrelativistic quantum critical systems at finite density
Fitzpatrick:2014efa ; Fitzpatrick:2014cfa . Notably, the example of fermions
near a Fermi surface interacting with a gapless boson studied therein also
exhibits UV/IR mixing.

$\displaystyle\frac{d\alpha_{0}(b)}{d\log(b)}\biggr{|}_{b=1}=\biggl{[}2-\frac{1}{\kappa}\log\biggl(\frac{2}{\tilde{\Lambda}}\biggr{missing})\biggr{]}\alpha_{0}~{}~{}.$
(52)

The key here is that, since $\tilde{\Lambda}\ll 1$, it follows that
$\alpha_{0}$ decays very quickly under the RG flow towards the infrared. Thus,
around the fractonic fixed point, this operator is irrelevant. This operator
does not destabilize the fractonic phase.

##  Appendix C Renormalization of the cosine gradient term

For the cosine gradient term, when $I=x$ and $f_{x}=a_{x}\partial_{x}h$, we
have

$\displaystyle g_{H}^{(xx)}(0)=\lim_{x^{\prime}\rightarrow
x}(2\pi)^{2}\bigl{\langle}a_{x}\partial_{x}h(x)\,a_{x}\partial_{x^{\prime}}h(x^{\prime})\bigr{\rangle}=(2\pi)^{2}\int\frac{d^{3}\mathbf{p}}{(2\pi)^{3}}\frac{(a_{x}p)^{2}}{\kappa\bigl{[}(pq)^{2}+k^{2}\bigr{]}}~{}~{}.$
(53)

Passing to dimensionless variables,

$\displaystyle g_{H}^{(xx)}(0)=\frac{\pi}{2\kappa}\int
d^{3}\tilde{\mathbf{p}}\,\frac{\tilde{p}^{2}}{(\tilde{p}\tilde{q})^{2}+\tilde{k}^{2}}=\frac{4\pi}{\kappa}\int_{>0}\frac{d\tilde{k}\,d\tilde{\ell}}{\tilde{k}^{2}+\tilde{\ell}^{2}}\int_{\tilde{\ell}}^{1}d\tilde{p}\,\tilde{p}=\frac{2\pi}{\kappa}\int_{>0}\frac{d\tilde{k}\,d\tilde{\ell}}{\tilde{k}^{2}+\tilde{\ell}^{2}}\,(1-\tilde{\ell}^{2})~{}~{}.$
(54)

By $\tilde{k}\leftrightarrow\tilde{\ell}$ symmetry, we can replace the factor
$1-\tilde{\ell}^{2}$ in the numerator with

$1-\tilde{\ell}^{2}\rightarrow
1-\frac{\tilde{k}^{2}+\tilde{\ell}^{2}}{2}~{}~{},$ (55)

and therefore,

$\displaystyle
g_{H}^{(xx)}(0)=\frac{2\pi}{\kappa}\int_{>0}d\tilde{k}\,d\tilde{\ell}\biggl{(}\frac{1}{\tilde{k}^{2}+\tilde{\ell}^{2}}-\frac{1}{2}\biggr{)}=\frac{\pi^{2}}{\kappa}\int_{\tilde{\Lambda}/b}^{\tilde{\Lambda}}d\tilde{\lambda}\biggl{(}\frac{1}{\tilde{\lambda}}-\frac{\tilde{\lambda}}{2}\biggr{)}=\frac{\pi^{2}}{\kappa}\biggl{[}\log(b)-\frac{\tilde{\Lambda}^{2}}{2}\biggl{(}1-\frac{1}{b^{2}}\biggr{)}\biggr{]}~{}~{}.$
(56)

The renormalized fugacity is

$\alpha_{x}(b)=\alpha_{x}\,b^{2-\frac{\pi^{2}}{2\kappa}}\,e^{-\frac{\pi^{2}}{4\kappa}\tilde{\Lambda}^{2}\bigl{(}1-\frac{1}{b^{2}}\bigr{)}}~{}~{},$
(57)

and the logarithmic derivative at $b=1$ evaluates to

$\displaystyle\frac{d\alpha_{x}(b)}{d\log(b)}\biggr{|}_{b=1}=\biggl{[}2-\frac{\pi^{2}}{2\kappa}\,(1+\tilde{\Lambda}^{2})\biggr{]}\alpha_{x}~{}~{}.$
(58)

We can drop the $\tilde{\Lambda}^{2}$ term as it is small by assumption:

$\displaystyle\frac{d\alpha_{x}(b)}{d\log(b)}\biggr{|}_{b=1}\approx\biggl{(}2-\frac{\pi^{2}}{2\kappa}\biggr{)}\alpha_{x}~{}~{}.$
(59)

Hence, there is a critical value

$\kappa_{c}=\frac{\pi^{2}}{4}~{}~{},$ (60)

at which $\alpha_{x}$ is marginal at this order, below which $\alpha_{x}$ is
irrelevant, and above which $\alpha_{x}$ becomes relevant. In other words,
this operator _does_ destabilize the fractonic phase when $\kappa>\kappa_{c}$.
This is precisely the process of dipole proliferation.

By reflection symmetry across the $x=y$ plane, the beta function for
$\alpha_{y}$ is the same as for $\alpha_{x}$, but with $x\rightarrow y$.

##  Appendix D Taking note of the discreteness of the gradients set by the
lattice scale

In the main text it was assumed that the lattice scale was small so that one
can turn lattice difference $\Delta_{i}$ into

$\displaystyle\Delta_{i}\rightarrow a_{i}\partial_{i}~{}~{},$ (61)

where $\Delta_{i}$ works on a function $j(\mathbf{x})$ as

$\displaystyle\Delta_{i}j(\mathbf{x})=j(\mathbf{x}+\mathbf{e}_{i})-j(\mathbf{x})~{}~{},$
(62)

where $\mathbf{e}_{i}$ is the unit vector in the $i$-direction. In Eq. (61)
repeated indices does not mean a contraction of indices. In Fourier space Eq.
(61) means

$\displaystyle\frac{2}{a_{i}}\sin(a_{i}k_{i}/2)\rightarrow k_{i}~{}~{}.$ (63)

We will now not use these simplifications and instead consider a discrete
model given by

$\displaystyle\begin{split}\mathcal{Z}&=\left[\prod_{\mathbf{x}}\int_{-\infty}^{\infty}dh\,\right]\text{exp}\Bigg{\\{}-\sum_{\mathbf{x}}\bigg{[}\frac{\kappa}{2}\biggl{(}\frac{\Delta_{x}\Delta_{y}}{a_{x}a_{y}}h\biggr{)}^{2}+\frac{\kappa}{2}\biggl{(}\frac{\Delta_{z}}{a_{z}}h\biggr{)}^{2}++2\sum_{I}\alpha_{I}\cos(2\pi
f_{I})\bigg{]}\Bigg{\\}}\,\,,\,\,f_{I}=\\{h,\Delta_{x}h,\Delta_{y}h\\}~{}~{},\end{split}$
(64)

For completeness we also discretized the $z$-direction in Eq. (64) but for
what follows we will take the $z$-direction to be continuous as we never come
close to the momentum cutoff $k^{max}=\pi/a_{z}$ so there is no reason to be
concerned that the approximation of Eq. (61) is problematic for this
direction. From these considerations it follows that the momentum shell one
should work with is given by

$\displaystyle\Lambda^{2}=\frac{16\sin^{2}(a_{x}p/2)\sin^{2}(a_{y}q/2)}{a^{2}_{x}a^{2}_{y}}+k^{2}~{}~{}.$
(65)

Furthermore, the simple Gaussian correlator in momentum space, which one needs
for the renormalization group analysis, is now modified by

$\displaystyle\langle
f_{0}(\mathbf{p})f_{0}(0)\rangle=\frac{1}{\kappa(p^{2}q^{2}+k^{2})}\rightarrow\frac{1}{\kappa}\frac{1}{\frac{16\sin^{2}(a_{x}p/2)\sin^{2}(a_{y}q/2)}{a^{2}_{x}a^{2}_{y}}+k^{2}}~{}~{}.$
(66)

Similarly

$\displaystyle\langle f_{x}(\mathbf{p})f_{x}(0)\rangle$
$\displaystyle=\frac{a_{x}^{2}p^{2}}{\kappa(p^{2}q^{2}+k^{2})}\rightarrow\frac{1}{\kappa}\frac{4\sin^{2}(a_{x}p/2)}{\frac{16\sin^{2}(a_{x}p/2)\sin^{2}(a_{y}q/2)}{a^{2}_{x}a^{2}_{y}}+k^{2}}~{}~{},$
(67a) $\displaystyle\langle f_{y}(\mathbf{p})f_{y}(0)\rangle$
$\displaystyle=\frac{a_{y}^{2}q^{2}}{\kappa(p^{2}q^{2}+k^{2})}\rightarrow\frac{1}{\kappa}\frac{4\sin^{2}(a_{y}q/2)}{\frac{16\sin^{2}(a_{x}p/2)\sin^{2}(a_{y}q/2)}{a^{2}_{x}a^{2}_{y}}+k^{2}}~{}~{}.$
(67b)

We are most interested in finding the quantitative effect of the
discretization on the critical value for the coefficient $\kappa$ above which
a transition takes place. For this we use Eq. (67a) to find

$\displaystyle g^{(xx)}_{H}(0)$
$\displaystyle=\frac{(2\pi)^{2}}{\kappa}\int_{\Lambda\rightarrow\Lambda/b}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}}\frac{4\sin^{2}(a_{x}p/2)}{\frac{16\sin^{2}(a_{x}p/2)\sin^{2}(a_{y}q/2)}{a^{2}_{x}a^{2}_{y}}+k^{2}}~{}~{}.$
(68)

Again, we pass to dimensionless variables as in Eq. (43) so that Eq. (68)
turns into

$\displaystyle g^{(xx)}_{H}(0)$
$\displaystyle=\frac{\pi^{3}}{2\kappa}\int_{\Lambda\rightarrow\Lambda/b}d^{3}\mathbf{\tilde{p}}\frac{4\sin^{2}(\pi\tilde{p}/2)}{16\sin^{2}(\pi\tilde{p}/2)\sin^{2}(\pi\tilde{q}/2)+\pi^{4}\tilde{k}^{2}}~{}~{}.$
(69)

It is difficult to evaluate Eq. (69) analytically because we have to take note
of the momentum shell given in Eq. (65). We therefore use two simplifications,
whose validity is justified by noting that upon reversing Eq. (61), one
reproduces the result found in the main text, which is free of these
approximations. Firstly, we consider a momentum shell along the $pq$-plane, so
that we can immediately integrate out $k$. This means the dimensionless shells
are given by

$\displaystyle\tilde{\Lambda}=16\sin^{2}(\pi\tilde{p}/2)\sin^{2}(\pi\tilde{q}/2)~{}~{}.$
(70)

When $\tilde{\Lambda}\ll 1$, the shell of a single quadrant can be
approximated by two orthogonally crossing (see fig. 2) ethanlake .

Figure 2: In Eq. (71) we simplify the momentum shell on the the upper right
quadrant of the $pq$-plane by turning it into two orthogonally crossing
strips. In both cases the width of the shell is $\Omega\varepsilon$, with
$\varepsilon=b-1$ and $\Omega$ being a positive constant which is unimportant
as it will drop out in the integration.

For either strip, we can linearize one sine, as for this strip either $p$ or
$q$ is always small. We will only be concerned with the strip where
$\tilde{q}$ is small, as the part of the integrand where $\tilde{p}$ is small
is negligible. Symmetry allows one to consider only this one quadrant and
multiply by $4$ to obtain the full integral. We can then integrate $\tilde{q}$
over the width of the strip to find

$\displaystyle\begin{split}g^{(xx)}_{H}(0)&\approx\frac{4\pi}{\kappa}\int_{0}^{1}d\tilde{p}\sin(\pi\tilde{p}/2)\log(b)=\frac{8}{\kappa}\log(b)~{}~{}.\end{split}$
(71)

In Eq. (71) we put the lower integration bound for $\tilde{p}$ to zero, as the
integrand is such that for $\tilde{\Lambda}\ll 1$ the additional contribution
that one obtains from taking the accurate lower integration bound instead is
negligible. From Eqs. (71) it follows that the flow of $\alpha_{x}$ is given
by

$\displaystyle\frac{\partial\alpha_{x}}{\partial\log(b)}=2\biggl{(}1-\frac{2}{\kappa}\biggr{)}\alpha_{x}~{}~{}.$
(72)

So one learns that the critical point for this discretized computation is
given by

$\displaystyle\kappa^{(L)}_{c}=2~{}~{}.$ (73)

## References

  * (1) Z. Zhou, X.-F. Zhang, F. Pollmann, and Y. You, “Fractal Quantum Phase Transitions: Critical Phenomena Beyond Renormalization,” arXiv:2105.05851 [cond-mat.str-el].
  * (2) Y. You and R. Moessner, “Fractonic plaquette-dimer liquid beyond renormalization,” arXiv:2106.07664 [cond-mat.str-el].
  * (3) E. Lake, “Renormalization group and stability in the exciton Bose liquid,” Physical Review B 105 no. 7, (Feb., 2022) 075115. https://link.aps.org/doi/10.1103/PhysRevB.105.075115.
  * (4) J. M. Kosterlitz and D. J. Thouless, “Ordering, metastability and phase transitions in two-dimensional systems,” Journal of Physics C: Solid State Physics 6 no. 7, (Apr, 1973) 1181–1203. https://doi.org/10.1088/0022-3719/6/7/010.
  * (5) V. Berezinskii, “Destruction of long-range order in one-dimensional and two-dimensional systems having a continuous symmetry group I. Classical systems,” Sov. Phys. JETP 32 no. 3, (1971) 493–500.
  * (6) P. Kasteleyn, “The statistics of dimers on a lattice: I. The number of dimer arrangements on a quadratic lattice,” Physica 27 no. 12, (1961) 1209–1225. https://www.sciencedirect.com/science/article/pii/0031891461900635.
  * (7) P. W. Kasteleyn, “Dimer Statistics and Phase Transitions,” Journal of Mathematical Physics 4 no. 2, (Feb, 1963) 287–293. https://doi.org/10.1063/1.1703953.
  * (8) M. E. Fisher and J. Stephenson, “Statistical Mechanics of Dimers on a Plane Lattice. II. Dimer Correlations and Monomers,” Phys. Rev. 132 (Nov, 1963) 1411–1431. https://link.aps.org/doi/10.1103/PhysRev.132.1411.
  * (9) D. S. Rokhsar and S. A. Kivelson, “Superconductivity and the Quantum Hard-Core Dimer Gas,” Phys. Rev. Lett. 61 (Nov, 1988) 2376–2379. https://link.aps.org/doi/10.1103/PhysRevLett.61.2376.
  * (10) C. L. Henley, “Relaxation time for a dimer covering with height representation,” Journal of Statistical Physics 89 no. 3, (Nov, 1997) 483–507. https://doi.org/10.1007/BF02765532.
  * (11) C. L. Henley, “From classical to quantum dynamics at Rokhsar–Kivelson points,” Journal of Physics: Condensed Matter 16 no. 11, (Mar, 2004) S891–S898. https://doi.org/10.1088/0953-8984/16/11/045.
  * (12) R. Moessner and K. S. Raman, “Quantum Dimer Models,” in Introduction to Frustrated Magnetism: Materials, Experiments, Theory, C. Lacroix, P. Mendels, and F. Mila, eds., pp. 437–479.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2011.  https://doi.org/10.1007/978-3-642-10589-0_17.
  * (13) R. Moessner and S. L. Sondhi, “Resonating Valence Bond Phase in the Triangular Lattice Quantum Dimer Model,” Phys. Rev. Lett. 86 (Feb, 2001) 1881–1884. https://link.aps.org/doi/10.1103/PhysRevLett.86.1881.
  * (14) C. Castelnovo, R. Moessner, and S. L. Sondhi, “Magnetic monopoles in spin ice,” Nature 451 no. 7174, (Jan, 2008) 42–45. https://doi.org/10.1038/nature06433.
  * (15) C. S. O. Yokoi, J. F. Nagle, and S. R. Salinas, “Dimer pair correlations on the brick lattice,” Journal of Statistical Physics 44 no. 5, (Sep, 1986) 729–747. https://doi.org/10.1007/BF01011905.
  * (16) M. Pretko, “Generalized electromagnetism of subdimensional particles: A spin liquid story,” Physical Review B 96 no. 3, (July, 2017) 035119. http://link.aps.org/doi/10.1103/PhysRevB.96.035119.
  * (17) M. Pretko, “Subdimensional particle structure of higher rank U(1) spin liquids,” Physical Review B 95 no. 11, (Mar., 2017) 115139. https://link.aps.org/doi/10.1103/PhysRevB.95.115139.
  * (18) N. Seiberg and S.-H. Shao, “Exotic $U(1)$ symmetries, duality, and fractons in 3+1-dimensional quantum field theory,” SciPost Physics 9 no. 4, (Oct., 2020) 046. https://scipost.org/10.21468/SciPostPhys.9.4.046.
  * (19) I. Herbut, A Modern Approach to Critical Phenomena.  Cambridge University Press, 2007.
  * (20) D. Tong, “Lectures on Statistical Field Theory,” 2017.
  * (21) Y. You, Z. Bi, and M. Pretko, “Emergent fractons and algebraic quantum liquid from plaquette melting transitions,” Phys. Rev. Research 2 (Feb, 2020) 013162. https://link.aps.org/doi/10.1103/PhysRevResearch.2.013162.
  * (22) 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 no. 10, (Oct, 2021) 102301. https://doi.org/10.1063%2F5.0060808.
  * (23) Y. You, J. Bibo, F. Pollmann, and T. L. Hughes, “Fracton Critical Point in Higher-Order Topological Phase Transition,” arXiv:2008.01746 [cond-mat.str-el].
  * (24) 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?,” arXiv:2101.01724 [cond-mat.str-el].
  * (25) J. V. José, L. P. Kadanoff, S. Kirkpatrick, and D. R. Nelson, “Renormalization, vortices, and symmetry-breaking perturbations in the two-dimensional planar model,” Phys. Rev. B 16 (Aug, 1977) 1217–1241. https://link.aps.org/doi/10.1103/PhysRevB.16.1217.
  * (26) A. Paramekanti, L. Balents, and M. P. A. Fisher, “Ring exchange, the exciton Bose liquid, and bosonization in two dimensions,” Phys. Rev. B 66 (Aug, 2002) 054526. https://link.aps.org/doi/10.1103/PhysRevB.66.054526.
  * (27) C. Xu and M. P. A. Fisher, “Bond algebraic liquid phase in strongly correlated multiflavor cold atom systems,” Phys. Rev. B 75 (Mar, 2007) 104428. https://link.aps.org/doi/10.1103/PhysRevB.75.104428.
  * (28) I. F. Herbut, B. H. Seradjeh, S. Sachdev, and G. Murthy, “Absence of U(1) spin liquids in two dimensions,” Physical Review B 68 no. 19, (Nov, 2003) . http://dx.doi.org/10.1103/PhysRevB.68.195110.
  * (29) I. F. Herbut and B. H. Seradjeh, “Permanent Confinement in Compact QED3 with Fermionic Matter,” Physical Review Letters 91 no. 17, (Oct, 2003) . https://doi.org/10.1103%2Fphysrevlett.91.171601.
  * (30) Z. Zhai and L. Radzihovsky, “Two-dimensional melting via sine-Gordon duality,” Physical Review B 100 no. 9, (Sep, 2019) . http://dx.doi.org/10.1103/PhysRevB.100.094105.
  * (31) J. M. Kosterlitz and D. J. Thouless, “Long range order and metastability in two dimensional solids and superfluids. (Application of dislocation theory),” Journal of Physics C: Solid State Physics 5 no. 11, (Jun, 1972) L124–L126. https://doi.org/10.1088/0022-3719/5/11/002.
  * (32) B. I. Halperin and D. R. Nelson, “Theory of Two-Dimensional Melting,” Phys. Rev. Lett. 41 (Jul, 1978) 121–124. https://link.aps.org/doi/10.1103/PhysRevLett.41.121.
  * (33) A. P. Young, “Melting and the vector Coulomb gas in two dimensions,” Phys. Rev. B 19 (Feb, 1979) 1855–1866. https://link.aps.org/doi/10.1103/PhysRevB.19.1855.
  * (34) A. L. Fitzpatrick, S. Kachru, J. Kaplan, S. Raghu, G. Torroba, and H. Wang, “Enhanced Pairing of Quantum Critical Metals Near d=3+1,” Phys. Rev. B 92 no. 4, (2015) 045118, arXiv:1410.6814 [cond-mat.str-el].
  * (35) A. L. Fitzpatrick, G. Torroba, and H. Wang, “Aspects of Renormalization in Finite Density Field Theory,” Phys. Rev. B 91 (2015) 195135, arXiv:1410.6811 [cond-mat.str-el].


* [27] Benjamin Assarf, Ewgenij Gawrilow, Katrin Herr, Michael Joswig, Benjamin Lorenz, Andreas Paffenholz, and Thomas Rehn.  Computing convex hulls and counting integer points with polymake, 2015.
  * [28] Ewgenij Gawrilow and Michael Joswig.  polymake: a Framework for Analyzing Convex Polytopes.  2000\.
  * [29] Klaus Hepp.  Proof of the Bogoliubov-Parasiuk theorem on renormalization.  Comm. Math. Phys., 2(1):301–326, Dec 1966.
  * [30] W. Zimmermann.  Convergence of Bogoliubov’s method of renormalization in momentum space.  Comm. Math. Phys., 15(3):208–234, Sep 1969.
  * [31] K.G. Chetyrkin and F.V. Tkachov.  Infrared r-operation and ultraviolet counterterms in the ms-scheme.  Physics Letters B, 114(5):340–344, 1982.
  * [32] K.G. Chetyrkin and V.A. Smirnov.  R-operation corrected.  Physics Letters B, 144(5):419–424, 1984.
  * [33] Franz Herzog and Ben Ruijl.  The R∗-operation for Feynman graphs with generic numerators.  JHEP, 05:037, 2017.
  * [34] Erik Panzer.  Algorithms for the symbolic integration of hyperlogarithms with applications to feynman integrals.  Computer Physics Communications, 188:148–166, March 2015.
  * [35] Máté L. Telek.  Geometry of the signed support of a multivariate polynomial and descartes’ rule of signs, 2023.
  * [36] Eduardo Cattani, Alicia Dickenstein, and Bernd Sturmfels.  Residues and resultants.  arXiv preprint alg-geom/9702001, 1997.
  * [37] Johannes M. Henn.  What can we learn about qcd and collider physics from n=4 super yang–mills?  Annual Review of Nuclear and Particle Science, 71(1):87–112, September 2021.
  * [38] Johannes M. Henn, Martín Lagares, and Shun-Qing Zhang.  Integrated negative geometries in abjm.  Journal of High Energy Physics, 2023(5), May 2023.
  * [39] Nima Arkani-Hamed, Johannes Henn, and Jaroslav Trnka.  Nonperturbative negative geometries: amplitudes at strong coupling and the amplituhedron.  Journal of High Energy Physics, 2022(3), March 2022.
  * [40] N. Arkani-Hamed, J. Bourjaily, F. Cachazo, and J. Trnka.  Local integrals for planar scattering amplitudes.  Journal of High Energy Physics, 2012(6), June 2012.
  * [41] Giulio Gambuti, David A. Kosower, Pavel P. Novichkov, and Lorenzo Tancredi.  Finite feynman integrals, 2023.
  * [42] Francis C. S. Brown.  On the periods of some feynman integrals, 2010.
  * [43] Francis Brown.  The Massless higher-loop two-point function.  Commun. Math. Phys., 287:925–958, 2009.
  * [44] Michael Borinsky.  Tropical Monte Carlo quadrature for Feynman integrals.  Annales de l’Institut Henri Poincaré D, jan 2023.
  * [45] Michael Borinsky, Henrik J. Munch, and Felix Tellander.  Tropical Feynman integration in the Minkowski regime.  Computer Physics Communications, 292:108874, nov 2023.
  * [46] Nima Arkani-Hamed, Hadleigh Frost, Giulio Salvatori, Pierre-Guy Plamondon, and Hugh Thomas.  All loop scattering for all multiplicity, 2023.
  * [47] Nima Arkani-Hamed, Carolina Figueiredo, Hadleigh Frost, and Giulio Salvatori.  Tropical Amplitudes For Colored Lagrangians.  2 2024.
  * [48] Nima Arkani-Hamed, Qu Cao, Jin Dong, Carolina Figueiredo, and Song He.  Scalar-Scaffolded Gluons and the Combinatorial Origins of Yang-Mills Theory.  12 2023.
  * [49] Nima Arkani-Hamed, Qu Cao, Jin Dong, Carolina Figueiredo, and Song He.  NLSM $\subset$ Tr$(\phi^{3})$.  1 2024.
  * [50] Alok Laddha and Amit Suthar.  Positive Geometries, Corolla Polynomial and Gauge Theory Amplitudes.  5 2024.
  * [51] Shounak De, Andrzej Pokraka, Marcos Skowronek, Marcus Spradlin, and Anastasia Volovich.  Surfaceology for Colored Yukawa Theory.  6 2024.