# RG and Stability in the Exciton Bose Liquid

Ethan Lake<EMAIL_ADDRESS>Department of Physics, Massachusetts Institute of
Technology, Cambridge, MA, 02139

###### Abstract

The exciton Bose liquid (EBL) is a hypothesized phase of bosons in 2+1D which
possesses a dispersion that is gapless along the coordinate axes in momentum
space. The low energy theory of the EBL involves modes on all length scales,
extending all the way down to the lattice spacing. In this note, we discuss an
RG scheme that can be used to address the stability of this and related phases
of matter. We find that in the absence of an extensively large symmetry group,
realizing the simplest formulation of the EBL always requires fine-tuning.
However, we also argue that the addition of certain marginal interactions can
be used to realize a stable phase, without the need for fine-tuning. A simple
generalization to 3+1D is also discussed.

##  I Introduction

Most fixed points of renormalization group (RG) flows are characterized by the
absence of any intrinsic length scale. At such fixed points the low energy
physics is scale-invariant (and often conformally invariant), with all the
universal data of the fixed point being determined by pure dimensionless
numbers.

One important counterexample is a Fermi liquid. Fermi liquids can be realized
as gapless endpoints of RG flows, but they are not scale-invariant: the Fermi
momentum $k_{F}$ plays a crucial role in determining correlation functions,
and low energy degrees of freedom live at length scales all the way down to
$1/k_{F}$, which is usually comparable to the lattice spacing.

Another less well-known example of a fixed point with an intrinsic length
scale is the exciton Bose liquid (EBL), which was proposed in Paramekanti et
al. (2002) and which has recently arisen in a diverse array of different
contexts (e.g. You et al. (2020a, b); Gorantla et al. (2021a); Seiberg and
Shao (2020a)). The EBL is a phase of bosons in 2+1D possessing a dispersion of
the form

$\varepsilon_{\mathbf{k}}\propto|\sin(k_{x}a/2)\sin(k_{y}a/2)|,$ (1)

where $a$ is the lattice spacing. The most important aspect of (1) is that
$\varepsilon_{\mathbf{k}}=0$ along the coordinate axes in momentum space. This
means that the theory possesses low energy modes on all length scales down to
the lattice spacing, and consequently exhibits “UV-IR mixing” Gorantla et al.
(2021a); Seiberg and Shao (2020a). As such, any theory capturing the universal
physics of this phase of matter cannot be fully scale-invariant, as it must
know about the scale $a$.

The EBL was originally proposed Paramekanti et al. (2002) to be a stable phase
of matter, viz. one which does not require fine-tuning to be realized. This
claim was however later disputed Tay et al. (2011). In order to have a clear
answer to the question of whether or not the EBL is stable (in the presence of
a given group of microscopic symmetries), one needs to construct an RG scheme
that determines which perturbations affect the universal physics. As far as
the author is aware, there does not seem to be a discussion of how to do this
in the literature (at least when the thermodynamic limit is taken), and in
this note we will attempt to fill in this gap.

In order to address the question of stability, we need to understand how to
perform RG in a system where low energy modes live at all length scales. RG is
often described as a procedure involving coarse-graining in space, a la
Kadanoff’s original spin-blocking procedure Kadanoff (1966). There is however
no unique way of performing RG, and the most useful scheme will depend on
context.

In general, a useful RG scheme is one which eliminates non-universal degrees
of freedom, namely those which are not necessary for describing the low energy
physics of the system. These non-universal degrees of freedom may or may not
be associated with short distance scales. In particular, associating RG flow
with a successive elimination of short distance degrees of freedom is only
appropriate in problems where things happening at short distances also live at
high energies, which is not the case for the EBL. It is therefore misleading
to describe the EBL and related models as being “beyond renormalization” You
and Moessner (2021); You et al. (2021)111One would certainly not use these
words when discussing a Fermi liquid, for example.; rather, such models simply
mandate an RG scheme which does not proceed by eliminating short-distance
modes. The purpose of this note is to construct an appropriate RG scheme, and
determine whether or not the EBL is in fact stable.

An outline of the remainder of this note is as follows. In section II, we give
a brief review of the aspects of the EBL which will be relevant in the
following sections. Section III describes our approach to performing RG, and
in IV we use this approach to analyze the stability of the EBL. We will find
that the simplest version of the EBL can only be realized with fine-tuning,
but that a certain choice of marginal “Landau parameters” can be made such
that a stable phase seems likely to exist. In section V we briefly describe a
generalization to a related 3+1D model, and we conclude in section VI.

##  II The EBL fixed point

In this introductory section we recapitulate the physics of the EBL from a
perspective that will be useful in subsequent sections. In all of what follows
we will be working in a setting appropriate for doing condensed matter
physics: we will be in (continuous) imaginary time, on a spatial lattice with
a finite lattice spacing $a$, and will be performing all calculations in the
thermodynamic limit. In this limit the system size $L$ is sent to infinity,
while $a$ is held fixed. For more background information and a discussion of
other types of starting assumptions, see e.g. Paramekanti et al. (2002); Tay
et al. (2011); Seiberg and Shao (2020a); Gorantla et al. (2021a). Working in
the thermodynamic limit is particularly important; if we were to instead use a
continuum limit in which $a$ is sent to zero, our conclusions about stability
would change, in line with the discussion of Ref. Gorantla et al. (2021a). We
will have more detailed comments to make about this issue later on.

The EBL is a system of bosons on a 2+1D square lattice at average density
$\mkern 1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5mu$. In most of
what follows we will take $\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5mu$ to be some generic
(incommensurate) value. The dynamics of the bosons is assumed to be dominated
by an onsite repulsion and a ring exchange hopping term, with the most
important terms in the microscopic Hamiltonian schematically of the form

$\displaystyle H$
$\displaystyle\sim\mathcal{K}\sum_{i}b(\mathbf{r}_{i})b^{\dagger}(\mathbf{r}_{i}+a{\mathbf{\hat{x}}})b(\mathbf{r}_{i}+a{\mathbf{\hat{x}}}+a{\mathbf{\hat{y}}})b^{\dagger}(\mathbf{r}_{i}+a{\mathbf{\hat{y}}})+h.c.$
(2) $\displaystyle\qquad+\mathcal{U}\sum_{i}(n(\mathbf{r}_{i})-\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5mu)^{2},$

where $n(\mathbf{r}_{i})$ is the boson number operator, $\mathbf{r}_{i}$ runs
over lattice sites, and $a$ is the lattice spacing.

One important feature of the ring-exchange term is that it separately
conserves the number of bosons along every row and column of the lattice. In
the absence of any other boson hopping terms in the Hamiltonian such as222Here
and in the following we will be abusing notation by letting $a$ stand for both
the lattice spacing, and, when appearing in a sum, a coordinate index
$a\in\\{x,y\\}$.
$\sum_{a}b^{\dagger}(\mathbf{r}_{i})b(\mathbf{r}_{i}+a{\mathbf{\hat{a}}})$,
the Hamiltonian thus possesses a gigantic group of subsystem symmetries, with
the boson number along each row and column of the lattice being separately
conserved. In what follows we will never include this subsystem symmetry as
part of our microscopic symmetry group, since a microscopic boson Hamiltonian
with this symmetry group requires a large amount of fine-tuning. Rather, we
will always imagine that the Hamiltonian above includes terms with small bare
coefficients which break the subsystem symmetry. Part of the task at hand is
to determine whether or not such terms are relevant (in the technical sense).
The actual microscopic symmetry group we will work with in this note will at
most consist of overall $U(1)$ boson number conservation, translation
symmetry, and the discrete symmetries of the square lattice. In fact, for the
purposes of the points we are trying to make, none of these symmetries are
essential, and we will eventually relax them in subsequent sections.

An analysis of the problem defined by the UV Hamiltonian (2) proceeds by using
a hydrodynamic description in terms of two compact fields $\phi,\theta$, which
keep track of fluctuations in the phase and density of the UV bosons,
respectively Paramekanti et al. (2002). This is done in a manner very similar
to what we would do when writing down a hydrodynamic description of
interacting bosons in 1+1D. The legitimacy of this approach can be justified a
posteriori by computing correlation functions using the $\phi,\theta$
description and noting their quasi 1+1D character, as well as by the fact that
the various ordered phases of the theory can be accessed from the hydrodynamic
description by turning on appropriate cosines of $\phi$ and $\theta$ (to be
discussed later).

In more detail, the hydrodynamic description works by performing a polar
decomposition on the microscopic boson annihilation operator by writing

$b=e^{i\phi}\sqrt{a^{2}\mkern 1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mu+\frac{1}{{2\pi}}\Delta_{x}\Delta_{y}\theta},$ (3)

where $\Delta_{a}$ denotes the dimensionless lattice gradient and $\phi$ is a
compact field keeping track of the boson phase. The subsystem symmetries
referred to above (which will always be broken microscopically) act as
$\phi(\mathbf{r})\mapsto\phi(\mathbf{r})+f(x)+g(y)$ for arbitrary functions
$f,g$. $\theta$ on the other hand is a field which keeps track of the
fluctuations in the boson density. It is defined on the sites of the dual
lattice, so that $\Delta_{x}\Delta_{y}\theta$ at a lattice site $\mathbf{r}$
can be written out as
$(\Delta_{x}\Delta_{y}\theta)(\mathbf{r})=\theta(\widetilde{\mathbf{r}})-\theta(\widetilde{\mathbf{r}}+a{\mathbf{\hat{x}}})+\theta(\widetilde{\mathbf{r}}+a{\mathbf{\hat{x}}}+a{\mathbf{\hat{y}}})-\theta(\mathbf{r}+a{\mathbf{\hat{y}}})$,
where $\widetilde{\mathbf{r}}=\mathbf{r}-(a/2,a/2)$.

As written in (3), $\theta$ must be constrained so that $a^{2}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mu+\frac{1}{{2\pi}}\Delta_{x}\Delta_{y}\theta$ has integer eigenvalues; as
usual this constraint will be enforced softly in the low-energy theory by
letting $\theta$ run over all real values, and adding cosines like
$\cos(q[{2\pi}\mkern 1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5muxy+\theta])$ to the low-energy action, with $q\in\mathbb{Z}$. The reason
for parametrizing the fluctuations in the density as
$\Delta_{x}\Delta_{y}\theta$ is because $\theta$ then determines the density
of quadrupolar ring-exchange configurations of bosons Paramekanti et al.
(2002), which given the form of the Hamiltonian are the most important density
fluctuations to keep track of. Note that this whole procedure is exactly
analogous to what we would do when studying interacting bosons in 1+1D, with
the only difference being that in the latter case we would replace (3) with
$b=e^{i\phi}\sqrt{a\mkern 1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mu+\frac{1}{{2\pi}}\Delta_{x}\theta}$.

As mentioned above, we will mostly be interested in scenarios where the
microscopic symmetries of interest are those of boson number conservation and
translation, together with $D_{4}$ point group symmetry. These act on the
hydrodynamic fields as Paramekanti et al. (2002)

$\displaystyle U(1)\,:$
$\displaystyle\phi(\mathbf{r})\mapsto\phi(\mathbf{r})+\lambda,\qquad\theta(\mathbf{r})\mapsto\theta(\mathbf{r})$
(4) $\displaystyle T_{\bm{\mu}}\,:$
$\displaystyle\phi(\mathbf{r})\mapsto\phi(\mathbf{r}+{\bm{\mu}}),$
$\displaystyle\theta(\mathbf{r})\mapsto\theta(\mathbf{r}+{\bm{\mu}})+{2\pi}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mu(\mu_{x}y+\mu_{y}x+\mu_{x}\mu_{y})$ $\displaystyle D_{4}\,:$
$\displaystyle\phi(\mathbf{r})\mapsto\phi(R_{\pi/4}^{n}R_{x}^{m}\mathbf{r}),$
$\displaystyle\theta(\mathbf{r})\mapsto(-1)^{n+m}\theta(R_{\pi/4}^{n}R_{x}^{m}\mathbf{r}),$

where $\lambda\in[0,2\pi)$, $T_{\bm{\mu}}$ is the operator which implements
translation through the vector ${\bm{\mu}}=(\mu_{x},\mu_{y})$,333The action of
$T_{\bm{\mu}}$ on $\theta$ can be understood by requiring that the boson
density $n=\mkern 1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mu+\frac{a^{2}}{{2\pi}}\Delta_{x}\Delta_{y}\theta$ transform as a density
under continuous translations, with $n(\mathbf{r})\mapsto
n(\mathbf{r}+{\bm{\mu}})(1+\nabla\cdot{\bm{\mu}})$ to first order in $\mu$. In
fact as written above this property only holds for transformations for which
$\nabla_{x}\nabla_{y}\mu^{y}=\nabla_{x}\nabla_{y}\mu^{x}=0$ (these
translations are more easily represented due to the form of the derivatives in
$\Delta_{x}\Delta_{y}\theta$), but restricting ourselves to these
transformations will be enough for the present purposes. Note also that the
transformation of $\theta$ is nonlinear, with the $\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5mu\mu^{x}\mu^{y}$ piece
ensuring that
$T_{{\bm{\mu}}}T_{{\bm{\nu}}}\theta=T_{{\bm{\mu}}+{\bm{\nu}}}\theta$. and we
have written a general element in $D_{4}$ in terms of a $\pi/4$ rotation
$R_{\pi/4}$ and a reflection about the $x$ axis $R_{x}$.

In terms of the hydrodynamic variables $\phi,\theta$, the appropriate action
for the Hamiltonian (2) can be written as

$\displaystyle S$ $\displaystyle=\int
d\tau\sum_{i}\,\Big{(}i\partial_{\tau}\phi\left(a^{2}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mu+\frac{1}{{2\pi}}\Delta_{x}\Delta_{y}\theta\right)+$ (5)
$\displaystyle\qquad\frac{R^{2}\varsigma}{{2\pi}}(1-\cos(\Delta_{x}\Delta_{y}\phi))+\frac{\varsigma}{{4\pi}R^{2}}(\Delta_{x}\Delta_{y}\theta)^{2}\Big{)}+\cdots,$

where the first term is the hydrodynamic representation of
$b^{\dagger}\partial_{\tau}b$, the sum over $i$ is over the implicit position
argument $\mathbf{r}_{i}$ of the $\phi,\theta$ fields, and where we have
chosen to parametrize the couplings in the Hamiltonian in terms of an energy
$\varsigma$ and a dimensionless constant $R^{2}$.444We will find the above
parametrization of the couplings in terms of $R^{2},\varsigma$ to be most
convenient in what follows. In terms of the standard notation used e.g. in
Paramekanti et al. (2002), we have $\varsigma\sim{2\pi}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mu\sqrt{\mathcal{U}\mathcal{K}},\,R^{2}=\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mu\sqrt{\mathcal{K}/\mathcal{U}}$. In (5) the terms in $\cdots$ include
other subdominant boson hopping terms such as
$b^{\dagger}(\mathbf{r}_{i}+a{\mathbf{\hat{a}}})b(\mathbf{r}_{i})+h.c.\xrightarrow{\sim}\cos(\Delta_{a}\phi)$,
as well as cosines like $\cos(q[{2\pi}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5muxy+\theta]),\cos(q[{2\pi}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5mux+\Delta_{y}\theta])$
etc. (with $q\in\mathbb{Z}$), which as mentioned above arise from softly
constraining the dimensionless boson number density to be integer-valued. The
terms written down explicitly in (5) preserve an additional subsystem symmetry
which acts on the $\theta$ fields (corresponding to the conservation of vortex
number on each row and column of the lattice), but this symmetry is broken
completely by the aforementioned cosines involving $\theta$.

As in the hydrodynamic analysis of interacting bosons in 1+1D, we can analyze
the low energy behavior of this system by first assuming that the cosines of
$\theta$ are small enough to allow $\theta$ to be integrated out via Gaussian
integration, producing an effective action in terms of $\phi$ alone. The
legitimacy of this step can then be determined a posteriori by using an RG
scheme to determine the relevance of the appropriate cosines. Doing this, we
then obtain555Here we are intentionally omitting the total derivative term
$\int d\tau\sum_{i}\mkern 1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mua^{2}\partial_{\tau}\phi$, which will not be important to keep track of
in what follows.

$S=\frac{R^{2}}{{2\pi}}\int
d\tau\sum_{i}\left(\frac{1}{2\varsigma}(\partial_{\tau}\phi)^{2}+\varsigma(1-\cos(\Delta_{x}\Delta_{y}\phi))\right)+\cdots,$
(6)

where the $\cdots$ again contain all the nonlinear interactions allowed by
symmetry and compactness of $\phi$.

Figure 1:  A color plot of the EBL dispersion, which is minimized along the
coordinate axes. The dashed black contour indicates a cutoff at the momentum
scale $\Lambda\ll 1/a$; the solid black contour indicates one at
$\Lambda-d\Lambda$. The modes lying between the two contours are integrated
out during an RG step.

Since we are assuming that the dynamics of the bosons is dominated by the
ring-exchange term, the modes which are relevant for describing the low energy
physics are those for which $\Delta_{x}\Delta_{y}\phi$ (but not necessarily
$\Delta_{a}\phi$) is small. This allows us to Taylor expand the cosine in (6)
to leading order, producing a quadratic action for $\phi$. Doing this, Fourier
transforming to momentum space, and then generalizing by letting
$R^{2},\varsigma$ become nontrivial functions of momentum, the quadratic part
of the action is then

$\displaystyle
S_{0}=\int_{\mathbf{k},\omega}\frac{R^{2}(\mathbf{k})}{{4\pi}\varsigma_{\mathbf{k}}}\left(\omega^{2}+\varepsilon_{\mathbf{k}}^{2}\right)|\phi(\omega,\mathbf{k})|^{2},$
(7)

where
$\int_{\mathbf{k},\omega}=\int\frac{d^{2}k}{({2\pi}/a)^{2}}\frac{d\omega}{{2\pi}}$,
and where the dispersion is

$\varepsilon_{\mathbf{k}}=4\varsigma_{\mathbf{k}}|\sin(k_{x}a/2)\sin(k_{y}a/2)|.$
(8)

This dispersion is shown in fig. 1, with its most salient feature being the
fact that it vanishes along the coordinate axes in momentum space.

Alternatively we may integrate out $\phi$ instead, producing a quadratic
action for $\theta$ that reads666The fact that $(\Delta_{x}\Delta_{y})^{2}$ is
not invertible doesn’t really matter, as the zero modes can be fixed away by
boundary conditions (just as in the case of inverting the Laplacian).

$\displaystyle\int_{\mathbf{k},\omega}\frac{1}{{4\pi}R^{2}(\mathbf{k})\varsigma_{\mathbf{k}}}(\omega^{2}+\varepsilon_{\mathbf{k}}^{2})|\theta(\omega,\mathbf{k})|^{2}.$
(9)

The duality between $\phi$ and $\theta$ is much the same as in the 1+1D case,
and simply sends $R^{2}(\mathbf{k})\leftrightarrow R^{-2}(\mathbf{k})$. This
free model, where all allowed cosines of $\theta$ are neglected (corresponding
to removing field configurations of $\phi$ containing vortices), is
essentially the same as the modified Villain XY-plaquette model introduced in
Gorantla et al. (2021b).

The functional form of $\varsigma_{\mathbf{k}}$ will turn out to not affect
any RG eigenvalues or other quantities of interest for the present discussion;
therefore in what follows we will simply assume that
$\varsigma_{\mathbf{k}}=\varsigma$ is independent of momentum.
$R^{2}(\mathbf{k})$ on the other hand does affect RG eigenvalues, and the
physics described by (7) depends in an essential way on its functional form.
In what follows we will allow $R^{2}(\mathbf{k})$ to be an arbitrary positive-
definite function that is smooth on momentum scales of order $1/a$, so that
the Fourier transform of $R^{2}(\mathbf{k})$ is local in real space.

##  III RG procedure

In this section we set up an RG scheme that can be used to address the
stability of the Gaussian action (7). The general approach we will take to RG
in the EBL is inspired by Shankar’s treatment of RG in Fermi liquids Shankar
(1994), and will essentially follow the procedure worked out in Lake et al.
(2021); Lake (to appear).

As mentioned in the introduction, the point of an RG analysis is to isolate
the universal physics, regardless of the length scales involved. In the
present setting, the universal physics of the fixed point is determined by the
modes living near the coordinate axes in momentum space; small modifications
to the dispersion in regions far away from the coordinate axes should
therefore be counted as irrelevant within any useful RG scheme. We therefore
impose a cutoff $\Lambda$ in momentum space by restricting to modes of $\phi$
and $\theta$ with momentum such that

$\varepsilon_{\mathbf{k}}\leq\varsigma\eta^{2},$ (10)

where we have defined the small parameter

$\eta\equiv a\Lambda\ll 1.$ (11)

An illustration of the cutoff imposed by (10) is shown in fig. 1. The fact
that $\eta\ll 1$ means that $\Delta_{x}\Delta_{y}\phi$ is always small, in
accordance with our assumption that we may get away with Taylor expanding the
cosine in (6).

The perturbations to the fixed point we will need to consider in our stability
analysis can all be expressed as777When $e^{i\mathcal{O}}$ has long-range
order, we will instead write the summand as
$1-\cos(\mathcal{O}(\mathbf{r}_{i},\tau))$.

$\delta S=g\eta^{d_{\mathcal{O}}}\varsigma\int
d\tau\,\sum_{i}\cos(\mathcal{O}(\mathbf{r}_{i},\tau)),$ (12)

where $\mathcal{O}$ is some equal-time polynomial in the $\phi,\theta$ fields,
$d_{\mathcal{O}}$ is a number whose determination will be discussed shortly,
and $g$ a small dimensionless coupling, of order $\eta^{0}$. Note that, as can
be shown using the scheme developed in this section, other terms involving
time derivatives such as $(\partial_{\tau}\Delta_{x}\phi)^{2}$ are either
marginal and go into modifying the function $R^{2}(\mathbf{k})$, or else are
irrelevant. Therefore in the following we will only consider terms without
time derivatives.

Even though the combination $g\varsigma\int d\tau$ is dimensionless, the
cutoff explicitly makes an appearance in (12) by way of the factor
$\eta^{d_{\mathcal{O}}}$. Determining the correct value to take for
$d_{\mathcal{O}}$ can be done by requiring that when evaluated on typical
field configurations, the integrand in (12) is of the same order in $\eta$ as
the kinetic term in the Gaussian fixed point action, viz. of order
$(\Delta_{x}\Delta_{y}\phi)^{2}\sim\eta^{4}$.888One can also determine
$d_{\mathcal{O}}$ by requiring that the perturbation $\delta S$ make
contributions to correlation functions / the free energy which are the same
order in $\eta$ as the appropriate quantities in the free theory. If
$d_{\mathcal{O}}$ is chosen to be larger than this value then $\delta S$ is
too small to have any affect as a perturbation, while if $d_{\mathcal{O}}$ is
chosen to be smaller than the effects of $\delta S$ are not perturbatively
small, which contradicts our assumptions. As we will see, $d_{\mathcal{O}}$
essentially determines the effective dimension that the “scaling” dimension of
$\mathcal{O}$ is to be compared to when determining the relevance of $\delta
S$ (scare quotes added here as we will not really be doing any scaling per
se—more on this shortly). Note that this effective dimension cannot simply be
determined by the number of spacetime dimensions in which $\cos(\mathcal{O})$
has power-law correlations—as was suggested in Paramekanti et al. (2002); Tay
et al. (2011)—since the nature of the problem means that our RG scheme cannot
involve any uniform rescaling of spacetime (real-space conformal perturbation
theory is rendered unworkably awkward for the same reason).

The correct value to take for $d_{\mathcal{O}}$ is determined on an operator-
by-operator basis. As an example, for typical field configurations in the
cutoff theory we have

$\Delta^{n}_{sx}\Delta^{m}_{py}\phi\sim\eta^{2\min(n,m)},$ (13)

where
$\Delta_{sa}\phi(\mathbf{r})=\phi(\mathbf{r})-\phi(\mathbf{r}+sa{\mathbf{\hat{a}}})$,
$s\in\mathbb{Z}$. Thus in order to have
$\eta^{d_{\Delta^{n}_{sx}\Delta^{m}_{py}\phi}}(1-\cos(\Delta^{n}_{sx}\Delta^{m}_{py}\phi))\sim(\Delta_{x}\Delta_{y}\phi)^{2}\sim\eta^{4}$
on typical field configurations, we must take

$d_{\Delta^{n}_{sx}\Delta^{m}_{py}\phi}=4(1-\min(n,m)).$ (14)

Therefore e.g. $d_{\Delta_{a}\phi}=4$, while an operator $\mathcal{O}$ which
is invariant under the row / column subsystem symmetries acting on $\phi$
(which can be written as (13) with both $n,m$ nonzero) has
$d_{\mathcal{O}}\leq 0$.

Now we move on to the determination of RG eigenvalues. In each RG step, we
first split up $\phi=\phi_{<}+\phi_{>}$ into fast and slow modes (and likewise
for $\theta$), with $\phi_{>}$ only containing modes satisfying

$(\eta^{\prime})^{2}<\varepsilon_{\mathbf{k}}/\varsigma<\eta^{2},$ (15)

and with $\phi_{<}$ containing the rest, where we have defined

$\eta^{\prime}\equiv\eta(1-dt)$ (16)

with $dt=d\Lambda/\Lambda\ll 1$ the RG time step. Note that the modes being
integrated out include modes of all frequencies — given the non-relativistic
nature of the problem and the lack of a need for a frequency cutoff when
calculating correlation functions, it is more natural to simply integrate out
all frequencies, thereby keeping the effective action local in time.

Integrating out the fast fields, to lowest order in $g$ the perturbation to
the slow field effective action is

$\delta S=g\varsigma\eta^{d_{\mathcal{O}}}\int
d\tau\sum_{i}\cos(\mathcal{O}_{<})e^{-\frac{1}{2}G_{\mathcal{O}_{>}}(0,{\mathbf{0}})}$
(17)

where $G_{\mathcal{O}_{>}}$ is the 2-point correlation function of
$\mathcal{O}_{>}$, with the decomposition
$\mathcal{O}=\mathcal{O}_{<}+\mathcal{O}_{>}$ induced from those of $\phi$ and
$\theta$. We then define the “scaling” dimension $\Delta_{\mathcal{O}}$ by

$G_{\mathcal{O}_{>}}(0,{\mathbf{0}})=4\Delta_{\mathcal{O}}dt+O(dt^{2}).$ (18)

The factor of $4$ here (a factor of 2 would be more normal) arises because
with this definition the power laws that appear in the correlation functions
of $\mathcal{O}$ are functions of spacetime distances to the power of
$2\Delta_{\mathcal{O}}$ (see appendix A). Rewriting the
$\eta^{d_{\mathcal{O}}}$ appearing in (17) in terms of $\eta^{\prime}$, we
then see that the mode integration effectively results in $g$ being replaced
by

$g^{\prime}=g(1+(d_{\mathcal{O}}-2\Delta_{\mathcal{O}})dt),$ (19)

so that the RG eigenvalue of $g$ is

$y_{g}=d_{\mathcal{O}}-2\Delta_{\mathcal{O}}.$ (20)

Therefore whether or not $g$ represents a relevant perturbation is determined
by comparing $2\Delta_{\mathcal{O}}$ to the effective dimension
$d_{\mathcal{O}}$.

Note that at no point have we rescaled coordinates so that the cutoff is
increased back to $\Lambda$; we have instead simply re-expressed $\eta$ in
terms of the new (reduced) cutoff. Due to the form of the dispersion a re-
scaling which returns the cutoff to its original value cannot be uniform in
momentum space, and as such must necessarily have a rather nasty
implementation in real space, which is where the hydrodynamic description of
the EBL is most naturally formulated (moreover, any such rescaling is
ultimately nothing more than a change of variables, and cannot by itself
contain any physical content).

##  IV Stability analysis

We now apply the general discussion of the preceding section to compute
scaling dimensions of operators in the EBL theory, with the goal of
determining the stability of the free fixed point (7).

###  IV.1 “Scaling” dimensions of operators

Before committing to a particular choice for the function $R^{2}(\mathbf{k})$,
let us make a few general comments. The operators added as perturbations to
the free fixed point that we are interested in can all be written as either
$\cos(\Phi_{q})$ or $\cos(\Theta_{q})$, where $\Phi_{q}$, $\Theta_{q}$
represent general integral linear combinations of $\phi,\theta$ fields,
respectively:

$\Phi_{q}=\sum_{i}q_{i}\phi(\mathbf{r}_{i}),\qquad\Theta_{q}=\sum_{i}q_{i}\theta(\mathbf{r}_{i}),\qquad
q_{i}\in\mathbb{Z}.$ (21)

The scaling dimensions of $\cos(\Phi_{q}),\ \cos(\Theta_{q})$ can be computed
as follows. First, consider $\cos(\Phi_{q})$ operators. Letting $dS_{\Lambda}$
denote the shell in momentum space containing fast modes with momenta
satisfying (15), the scaling dimension of $\cos(\Phi_{q})$ is extracted by
computing the fast-mode correlator as (working to leading order in $dt$)

$\displaystyle 4dt\Delta_{\Phi_{q}}$
$\displaystyle=\int_{dS_{\Lambda}}\frac{dk_{x}dk_{y}}{({2\pi}/a)^{2}}\int_{\mathbb{R}}d\omega\,\frac{|q(\mathbf{k})|^{2}}{R^{2}(\mathbf{k})}$
(22)
$\displaystyle\qquad\times\frac{1}{\omega^{2}/\varsigma+\varsigma(4|\sin(k_{x}a/2)\sin(k_{y}a/2)|)^{2}}$
$\displaystyle=\frac{1}{{4\pi}\Lambda^{2}}\int_{dS_{\Lambda}}dk_{x}\,dk_{y}\,\frac{|q(\mathbf{k})|^{2}}{R^{2}(\mathbf{k})},$

with $q(\mathbf{k})$ the lattice Fourier transform of $q_{i}$.

As mentioned above, $R^{2}(\mathbf{k})$ is assumed by spatial locality to be
the Fourier transform of a function which is localized on length scales below
$1/\Lambda^{2}a=1/\eta\Lambda$. This means that $R^{2}(\mathbf{k})$ can be
Taylor expanded about zero momentum when either $k_{x}$ or $k_{y}$ is much
less than $\pi/a$, in particular when $k_{x},k_{y}\sim\Lambda$. We furthermore
will only be interested in operators $\Phi_{q}$ which are themselves local on
scales below $1/\eta\Lambda$, so that $|q(\mathbf{k})|^{2}$ can be similarly
expanded. Now the integral over the shell $dS_{\Lambda}$ can be split into
regions where $k_{x}a\ll 1$ and $k_{y}$ ranges from $\Lambda$ to $\pi/a$ (for
which the shell $dS_{\Lambda}$ is defined via
$k_{x}=\Lambda^{2}a/2\sin(k_{y}a/2)$), and likewise for $k_{x}\leftrightarrow
k_{y}$. Performing appropriate Taylor expansions of $R^{2}(\mathbf{k})$ and
$q(\mathbf{k})$ in these regions, we may perform the integral over
$dS_{\Lambda}$ and write

$\displaystyle\Delta_{\Phi_{q}}$
$\displaystyle=\frac{a}{4\pi}\int_{\Lambda}^{\pi/a}\Big{(}\frac{dk_{x}\,|q(k_{x},\Lambda^{2}a/2\sin(k_{x}a/2))|^{2}}{R^{2}(k_{x},\Lambda^{2}a/2\sin(k_{x}a/2))\sin(ak_{x}/2)}$
(23)
$\displaystyle\qquad+\frac{dk_{y}\,|q(\Lambda^{2}a/2\sin(k_{y}a/2),k_{y})|^{2}}{R^{2}(\Lambda^{2}a/2\sin(k_{y}a/2),k_{y})\sin(ak_{y}/2)}\Big{)}$
$\displaystyle\approx\frac{a}{4\pi}\int_{\Lambda}^{\pi/a}\Big{(}\frac{dk_{x}\,|q(k_{x},0)|^{2}}{R^{2}(k_{x},0)\sin(ak_{x}/2)}$
$\displaystyle\qquad+\frac{dk_{y}\,|q(0,k_{y})|^{2}}{R^{2}(0,k_{y})\sin(ak_{y}/2)}\Big{)},$

where $\approx$ means equality up to terms suppressed by higher powers of
$\eta$.

There are several things to note about this expression. First, note that the
scaling dimension depends only on the values that $R^{2}(\mathbf{k})$ takes on
the coordinate axes. Therefore for the purposes of determining the stability
of the Gaussian fixed point we only need to know the function $R^{2}(k_{x},0)$
(which is equal to $R^{2}(0,k_{y})$ by the square lattice symmetry we have
assumed to be present for simplicity).

Secondly, note that $\Delta_{\Phi_{q}}$ diverges logarithmically as
$-\ln(\eta)$ unless $q({\mathbf{0}})=0$, since $R^{2}({\mathbf{0}})$ is a
assumed to be finite. This implies that $\Delta_{\Phi_{q}}$ diverges if
$\sum_{i}q_{i}\neq 0$, i.e. if $\Phi_{q}$ is charged under the global $U(1)$
boson number symmetry. Hence only operators which conserve total boson number
have a chance to be relevant.

Thirdly, note that the scaling dimension vanishes with $\eta\rightarrow 0$ if
$q(k_{x},0)=q(0,k_{y})=0$ for all $k_{x},k_{y}$. This condition is equivalent
to the condition that $\Phi_{q}$ be neutral under the row and column subsystem
symmetries which act on $\phi$. As was discussed near (13), any operator
neutral under both row and column subsystem symmetries must have
$d_{\mathcal{O}}\leq 0$. Together with the fact that such operators have
$\Delta_{\mathcal{O}}=O(\eta)\rightarrow 0$, this means that any perturbation
preserving the subsystem symmetries is guaranteed to be either irrelevant or
marginal.

Therefore any operator $\cos(\Phi_{q})$ which has a chance to be relevant must
both a) preserve the global $U(1)$ symmetry and b) break the subsystem
symmetries. From the discussion near (13) we see that such operators have
$d_{\mathcal{O}}=4$; as such their relevance is determined by comparing their
scaling dimensions to 2.

Above we have focused on $\cos(\Phi_{q})$ operators. The story for the
$\cos(\Theta_{q})$ operators is the same, with the only difference in the
calculation of $\Delta_{\Theta_{q}}$ being the replacement
$R^{2}(\mathbf{k})\rightarrow R^{-2}(\mathbf{k})$. In particular, any operator
with nonzero vortex number (such as $\cos(\theta)$) will be infinitely
irrelevant, while any operator invariant under the dual subsystem symmetries
(which count the vortex number in each row and column of the lattice) will be
either irrelevant or marginal.

From (23) is clear that the operators with the smallest nonzero scaling
dimensions must either have $q(k_{x},0)=0$ for all $k_{x}$ or $q(0,k_{y})=0$
for all $k_{y}$, but not both. Therefore when going about finding operators
which have the potential to destabilize the fixed point, we may restrict our
attention to operators which break one of the row / column subsystem
symmetries, but not both. Without loss of generality we may thus only consider
operators with nonzero $q(k_{x})\equiv q(k_{x},0)$. Combined with the fact
that the scaling dimensions of such operators depend only on $R^{2}(k_{x},0)$,
the calculation of the smallest scaling dimensions appearing in the operator
spectrum reduces to a one-dimensional optimization problem, considerably
simplifying the stability analysis.

When performing a numerical search for relevant operators, it is helpful to
further simplify things slightly. Since all $\Phi_{q}$ ($\Theta_{q}$)
operators we need to consider have zero (vortex) charge, implying that
$q(0)=0$, we may without loss of generality write $q(k)$ in terms of integers
$q_{l}$ defined on the links of a one-dimensional lattice as

$q(k)=\sum_{l\in\mathbb{Z}}e^{-ikla}q_{l}(1-e^{-ika}).$ (24)

In terms of the $q_{l}$, the scaling dimensions are

$\displaystyle\Delta_{\Phi_{q}}$
$\displaystyle=\frac{1}{\pi}\int_{0}^{\pi/a}dk\,\frac{\sin(ka/2)}{R^{2}(k_{x},0)}\sum_{l,l^{\prime}}q_{l}q_{l^{\prime}}\cos([l-l^{\prime}]ka)$
(25) $\displaystyle\Delta_{\Theta_{q}}$
$\displaystyle=\frac{1}{\pi}\int_{0}^{\pi/a}dk\,R^{2}(k_{x},0)\sin(ka/2)\sum_{l,l^{\prime}}q_{l}q_{l^{\prime}}\cos([l-l^{\prime}]ka).$

The Guassian fixed point parameterized by the given choice of
$R^{2}(\mathbf{k})$ is then stable provided that there are no nontrivial
choices of $\\{q_{l}\\}$ for which at least one of
$\Delta_{\Phi_{q}},\Delta_{\Theta_{q}}$ is less than $2$.

For now, we will only be interested in perturbations which respect translation
symmetry. From the translation action of (4), it is easy to check that at
generic incommensurate $\mkern 1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mu$, any translation-invariant operator $\cos(\Theta_{q})$ must have zero
vortex dipole moment999This statement applies to cosines of integral linear
combinations of $\theta$ fields. One may also use explicit coordinate
dependence to produce translation invariant cosines, e.g. as
$\cos({2\pi}\mkern 1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5muax-\Delta_{y}\theta)$. However as any finite-action field configuration
must asymptotically have $\Delta_{x}\Delta_{y}\theta\rightarrow 0$, such
cosines are always rapidly oscillating at large distances for incommensurate
$\mkern 1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5mu$, and as such
can be ignored Paramekanti et al. (2002). (such as e.g.
$\cos(\Delta_{a}^{2}\theta)$). Therefore translation invariance allows us to
restrict the $\\{q_{l}\\}$ in the second line of (25) to those integers
satisfying $\sum_{l}q_{l}=0$.

###  IV.2 Constant $R^{2}(\mathbf{k})$

Let us first consider the case where $R^{2}(k_{x},0)=R^{2}$ is a constant,
independent of momentum. In this case there turns out to always exist a
symmetry-allowed relevant perturbation.

To show this, we start by considering the simplest $\Phi_{q}$ operators which
preserve the total boson number, viz. exponentials of $\Delta_{x}\phi$ and
$\Delta_{y}\phi$. The dimension of these operators is, reading off from (25),

$\Delta_{\Delta_{x}\phi}=\frac{1}{\pi
R^{2}}\int_{0}^{\pi}dx\,\sin(x/2)=\frac{2}{\pi R^{2}}.$ (26)

As $d_{\Delta_{a}\phi}=4$ (so that $\Delta_{\Delta_{x}\phi}$ should be
compared to 2 when determining relevance), stability of the EBL fixed point
with constant $R^{2}$ requires

$R^{2}<1/\pi.$ (27)

On the other hand, consider translation-invariant operators $\Theta_{q}$ built
from the $\theta$ fields. The simplest such operators are those involving two
derivatives, which may be written as $\Delta_{sx}\Delta_{px}\theta$ for some
integers $s,p\neq 0$. These operators have scaling dimensions

$\Delta_{\Delta_{sx}\Delta_{px}\theta}=\frac{4R^{2}}{\pi}\int_{0}^{\pi}dx\,\frac{\sin^{2}(sx/2)\sin^{2}(px/2)}{\sin(x/2)}.$
(28)

The RHS is minimized when one of $s,p$ is equal to unity, with the other made
large. Specifically, if we set $p=1$ (such operators were previously
identified in Tay et al. (2011)), we find

$\Delta_{\Delta_{sx}\Delta_{x}\theta}=\frac{4R^{2}}{\pi(1-1/4s^{2})}.$ (29)

Strictly speaking, the most relevant operator is therefore
$\Delta_{(s\rightarrow\infty)x}\Delta_{x}\theta$. However, on general grounds
the bare coupling constant for
$\cos(\Delta_{(s\rightarrow\infty)x}\Delta_{x}\theta)$ will decay (usually
exponentially fast) with $s$. Since the variation of the scaling dimension
with $s$ is rather small, we will simply restrict our attention to $s=1$, with
the scaling dimension

$\Delta_{\Delta_{x}^{2}\theta}=\frac{16R^{2}}{3\pi}.$ (30)

This operator is therefore only irrelevant provided that

$R^{2}>\frac{3\pi}{8}.$ (31)

By examining (27) and (31), we see that no matter the value of $R^{2}$, there
is always a relevant perturbation. We now briefly discuss the nature of the RG
flows driven by these perturbations.

####  IV.2.1 $\cos(\Delta_{a}\phi)$ most relevant

If $\cos(\Delta_{a}\phi)$ is the most relevant perturbation and a term like
$\sum_{a}[1-\cos(\Delta_{a}\phi)]$ is added to the action, the flow drives the
model to a regime where $\Delta_{a}\phi\ll 1$ on typical field configurations,
so that the aforementioned cosine can be Taylor expanded. The resulting
$\sum_{a}(\Delta_{a}\phi)^{2}$ term lifts the degeneracy of the dispersion,
and in the IR we obtain the ordered phase of the 2+1D XY model. In particular,
the IR theory is massless, and the response to a background $U(1)$ gauge field
is that of a superconductor (unlike the response of the EBL phase described by
(7), which is insulating Paramekanti et al. (2002)).

####  IV.2.2 $\cos(\Delta^{2}_{a}\theta)$ most relevant

Now consider the case where $\cos(\Delta^{2}_{a}\theta)$ is the most relevant
perturbation. As was discussed above, technically speaking
$\cos(\Delta_{sa}\Delta_{a}\theta)$ is more relevant for larger $s$, but will
also appear with a bare coupling constant that is suppressed rather quickly
with large $s$. Since the effects of these operators for different small
values of $s$ are similar, we will simply assume that the most important
operator for determining the RG flow is $\cos(\Delta_{a}^{2}\theta)$. These
operators break down the group of subsystem symmetries acting on $\theta$ in
(9) to the subgroup corresponding to conservation of total vortex charge and
total vortex dipole moment (aka momentum).

In this case, when a term like $\sum_{a}[1-\cos(\Delta^{2}_{a}\theta)]$ is
added to the action, we flow to a regime where $\Delta_{a}\Delta_{b}\theta\ll
1$ for all $a,b=x,y$. After Taylor expanding, the free part of the action can
then be written in continuum notation as

$S_{0}=\int
d\tau\,d^{2}x\,\left((\partial_{\tau}\theta)^{2}+A\sum_{a}(\partial_{a}^{2}\theta)^{2}+B(\partial_{x}\partial_{y}\theta)^{2}\right)$
(32)

for constants $A,B$. This is essentially the quantum Lifshitz model with a
square lattice anisotropy,101010The usual quantum Lifshitz model on the square
lattice is generically unstable Vishwanath et al. (2004); Fradkin et al.
(2004) due to nonlinear terms involving derivatives like
$(\partial_{x}\theta)^{4}+(\partial_{y}\theta)^{4}$. In the present case such
terms are disallowed by translation invariance.  with the added proviso that,
by translation invariance, all terms must preserve the vortex dipole moment
(viz. must either involve time derivatives, or at least two spatial
derivatives). There is no longer any degeneracy of the dispersion along the
coordinate axes. However, since the terms in (32) are all quartic in spatial
derivatives, $e^{i\theta}$ has spatial correlation functions going
schematically like $\sim\exp\left(-\int
d^{2}k\,\frac{1-\cos(\mathbf{k}\cdot\mathbf{r})}{k^{2}}\right)$, and therefore
cannot develop long-range order. On the other hand, the operator
$e^{i\Delta_{a}\theta}$ does have long-range correlators, and consequently
develops long-range order, spontaneously breaking the conservation of vortex
dipole moment. Since $\Delta_{a}\theta$ is charged under translation, the
resulting state spontaneously breaks translation symmetry, and possesses a
gapless phonon mode. That this phase is gapless is in keeping with a general
analysis of the t’ Hooft anomalies present in the fixed point action (7), with
an anomaly being present unless either one subsystem symmetry is fully broken,
or both subsystem symmetries are broken down to their global $U(1)$ subgroups
Lam (private communication).

The exact type of charge ordering that occurs can be determined by writing
down a low-energy expression for the density operator. In the microscopic
boson model, the fluctuations in the number density are given by $\delta
n=\frac{1}{{2\pi}a^{2}}\Delta_{x}\Delta_{y}\theta$, as written down in (3). By
the time enough high-energy degrees of freedom are integrated so as to land us
in the low-energy hydrodynamic description, this expression will be
generically modified to include all other terms involving $\theta$ which are
allowed by symmetry. The constraints of $D_{4}$ symmetry and the particle-hole
symmetry occurring at half-filling mean that we may write (see e.g. Tay et al.
(2011); Xu and Fisher (2007) for related discussions)

$\displaystyle\delta n(\mathbf{r})$
$\displaystyle=\frac{1}{{2\pi}a^{2}}(\Delta_{x}\Delta_{y}\theta)(\mathbf{r})$
(33)
$\displaystyle+\sum_{q\in\mathbb{N}}\Upsilon_{q}\Delta_{x}\Delta_{y}\sin(q[\theta+{2\pi}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5muxy]))(\mathbf{r})$
$\displaystyle+\sum_{p\in\mathbb{N}}\Gamma_{p}\left((\Delta_{x}\sin(p[\Delta_{y}\theta+{2\pi}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5muax]))(\mathbf{r})+(x\leftrightarrow y)\right)$ $\displaystyle+\cdots,$

where the $\Upsilon_{q},\Gamma_{p}$ are non-universal coefficients, and the
$\cdots$ represent subleading contributions to $\delta n(\mathbf{r})$
involving $\theta$ fields living beyond the four dual lattice sites nearest to
$\mathbf{r}$. In the ordered phase we are interested in, we have $\langle
e^{i\Delta_{a}\theta}\rangle\neq 0$ and $\langle e^{i\theta}\rangle=0$. In
this case it is the components of the density involving $\Gamma_{p}$ which
order, and we may write

$\displaystyle\langle\delta n\rangle$ $\displaystyle\approx{2\pi}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mua^{2}\sum_{p\in\mathbb{N}}p\Gamma_{p}\big{(}\langle\cos(p\Delta_{x}\theta)\rangle\cos({2\pi}p\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5muay)$ (34)
$\displaystyle\qquad-\langle\sin(p\Delta_{x}\theta)\rangle\sin({2\pi}p\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5muay)\big{)}+(x\leftrightarrow y).$

###  IV.3 General $R^{2}(\mathbf{k})$

We now turn to asking whether or not there exist more complicated choices of
$R^{2}(k_{x},0)$ that yield phases stable against the condensation of either
particle or vortex dipoles.

In general, this problem can be formulated as a shortest vector problem, where
the task is to identify the length of the shortest vector in the integral
lattices definable from (25). This problem is generically very difficult
(especially as the dimension of the lattice increases; here the dimension goes
to infinity in the thermodynamic limit), and in the absence of more
sophisticated constructive approaches like those of Plamadeala et al. (2014),
we must resort to a brute-force numerical search. This is ameliorated somewhat
by the fact that a large number of choices for $\\{q_{l}\\}$ can be ruled out
from the beginning (for example, we may require $\gcd(\\{q_{l}\\})=1$ and
$\sum_{l}q_{l}\geq 0$ without loss of generality), but it still means that we
will only at best be able to provide suggestive evidence for (but not prove)
the existence of a stable phase.

In the following, we perform our numerical search by restricting ourselves to
operators with involve no more than $\sf{MaxBody}$ boson / vortex creation and
annihilation operators, with the given operators separated by no more than
$\sf{MaxRange}$ lattice points. To rigorously demonstrate stability, we would
need to make arguments for what happens as
$\mathsf{MaxBody},\mathsf{MaxRange}\rightarrow\infty$. Instead, we will simply
content ourselves with performing a stability analysis for a series of
increasing values of $\mathsf{MaxBody},\mathsf{MaxRange}$, and observing
whether or not any putative choices of $R^{2}(\mathbf{k})$ appear to be stable
as $\mathsf{MaxBody},\mathsf{MaxRange}$ are increased. Note that if the system
develops an instability only when $\mathsf{MaxBody}$ or $\mathsf{MaxRange}$
are increased past some large value, the EBL with the given choice of
$R^{2}(\mathbf{k})$ may be able to be regarded as stable in practice, as very
high body operators / those that act over a large number of sites will
generically enter the UV action with very small coefficients, preventing them
from appearing until the very latest stages of the RG flow (before which in
any real system the flow will be cut off either by finite size effects or by
nonzero temperature).

There are of course an infinite number of functions $R^{2}(k_{x},0)$ to try
when searching for a stable phase. One simple choice which seems to do a good
job is

$R(k_{x},0)=\lambda_{1}(1+\lambda_{2}\cos(ak_{x}))^{2},$ (35)

with $\lambda_{1}>0$ and $|\lambda_{2}|<1$. Note that this meets the
requirements of being positive and being the Fourier transform of a reasonably
localized function.

Figure 2:  Stability of the EBL for $R^{2}(k_{x},0)$ given by (35) and
$\mathsf{MaxBody}=6,\mathsf{MaxRange}=6$. Green regions indicate stability,
i.e. the absence of relevant operators. Left: regions where no relevant
$\cos(\Phi_{q})$ operators exist. Center: regions where no relevant
$\cos(\Theta_{q})$ operators exist. Right: regions where no relevant operator
of either type exists.

A plot showing the regions in $\lambda_{1},\lambda_{2}$ parameter space where
(35) yields a stable phase for $\mathsf{MaxBody},\mathsf{MaxRange}=6$ is shown
in 2. We see that there exist small regions of stability (green regions in the
rightmost panel of fig. 2), which are located near the regions where
$|\lambda_{2}|=1$. The fact that the regions of stability are located near the
border of the allowed parameter space is a common theme in these types of
problems Lake and Hermele (2021); Mukhopadhyay et al. (2001); Vishwanath and
Carpentier (2001).

Figure 3:  Size of the region of stability in Fig. 2 near
$-1<\lambda_{2}<-0.75$, $1<\lambda_{1}<4$, as a function of $\mathsf{MaxBody}$
and $\mathsf{MaxRange}$. Here “stability fraction” denotes the area in
$\lambda_{1}$-$\lambda_{2}$ parameter space of the stable region, relative to
the area when $\mathsf{MaxBody}=\mathsf{MaxRange}=4$.

Upon increasing $\mathsf{MaxBody},\mathsf{MaxRange}$, the region of stability
shrinks somewhat (especially at $\mathsf{MaxBody}=8$), but does not completely
disappear up to the largest values of $\mathsf{MaxBody},\mathsf{MaxRange}$ we
have numerically checked. The evolution of the size of the stable region with
$\mathsf{MaxBody},\mathsf{MaxRange}$ is shown in fig. 3. Whether or not these
curves should be viewed as extrapolating to a nonzero value in the
$\mathsf{MaxBody},\mathsf{MaxRange}\rightarrow\infty$ limit is a decision left
to the reader, but we remark again that even if the allowed region of
stability vanishes when $\mathsf{MaxBody}$ or $\mathsf{MaxRange}$ are very
large, the EBL in this parameter range may still be stable for all practical
purposes.

In the regions of instability, the exact RG flow will depend on the form of
the most relevant operators, as well as the strengths of their bare coupling
constants. In general though, we expect that the flow will be as in the
constant $R^{2}$ case, viz. either towards the superfluid phase of the 2+1D XY
model, or towards a translation-breaking crystalline state.

###  IV.4 Stability for other symmetry groups

So far we have restricted ourselves to perturbations which respect the
symmetries of translation and total boson number conservation (as well as
square lattice symmetry, though the latter is inessential). Boson number
conservation is inconsequential to our stability analysis, since as we have
seen operators which do not conserve boson number have infinite scaling
dimensions. However, the stability analysis will indeed change if we relax our
imposition of translation symmetry, or if we impose additional symmetries.

####  IV.4.1 No translation symmetry / commensurate density

Figure 4:  The same as in fig. 2, but working at integral densities /
including operators which break translation symmetry.

Suppose now that we ignore the requirement that perturbations to the fixed
point theory preserve translation symmetry. This then forces us to consider
perturbations involving any combinations of $\theta$ fields with zero vortex
number111111Operators with nonzero vortex number are still irrelevant, by the
reasoning given earlier. regardless of their dipole moment, such as e.g.
$\sum_{a}\cos(\Delta_{a}\theta)$. The same type of perturbations are allowed
if we keep translation symmetry but work at a commensurate density with
$a^{2}\mkern 1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mu\in\mathbb{Z}^{\geq 0}$, so that on average there are an integer number
of bosons per site. In this case the factors involving explicit coordinate
dependence which appear in cosines involving $\theta$ fields (such as the
${2\pi}\mkern 1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5mux$ in
$\cos(\Delta_{y}\theta+{2\pi}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5muax)$) can all be dropped
on account of ${2\pi}\mkern 1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5muax\in{2\pi}\mathbb{Z}$ for all lattice sites $x$. From the perspective of
stability, this case is therefore equivalent to the one where we ignore
translation symmetry.

In this setting, we may now consider arbitrary sets of integers $\\{q_{l}\\}$
in the second line of (25). Obviously any region of stability in this case
must be a proper subset of the region of stability found in the case where
translation symmetry was imposed microscopically. We find that for the choice
of $R^{2}(k_{x},0)$ in (35), the region of (apparent) stability is reduced but
not altogether eliminated, as shown in fig. 4.

Consider the regions of instability where the RG flow is not towards the
superfluid phase. If $a^{2}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5mu$ is not an integer, the
flow will be towards some state with a pattern of charge order determined by
the most relevant perturbation. If instead $a^{2}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5mu$ is integral, the flow
will be towards a translation-invariant Mott insulator with $a^{2}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5mu$ bosons per site. Either
way, the resulting phase will be massive, which is allowed by anomaly
constraints due to the dual vortex subsystem symmetry being completely broken
Lam (private communication).

If translation symmetry is imposed and $a^{2}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5mu$ is not an integer but
some relatively commensurate rational number, more cosines involving $\theta$
operators are allowed, even in the presence of translation symmetry. The
region of stability in this case will then be somewhere between the regions
shown in figs. 2 and 4, depending on the value of $a^{2}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern 1.5mu$. In the regions of
instability, if the most relevant operator is a cosine of $\theta$, the
resulting RG flow will generally be towards some sort of charge density wave;
see Paramekanti et al. (2002); Tay et al. (2011) for a detailed discussion.

####  IV.4.2 Dipole conservation

We may also consider a theory with a larger (but still finite-dimensional)
global symmetry group. One symmetry we may impose is global dipole
conservation, which maps $\phi\mapsto\phi+\alpha x+\beta y$ for constant
$\alpha,\beta$, and under which operators like $\cos(\Delta_{x}\phi)$ carry
charge. If we impose this symmetry in addition to translation, the region of
stability in fig. 2 can only increase.

Dipole conservation together with translation symmetry is however not enough
to render the theory with constant $R^{2}$ stable. Indeed, irrelevance of
$\cos(\Delta_{x}^{2}\theta)$ still requires that $R^{2}>3\pi/8$, while the
simplest dipole-neutral operator $\cos(\Delta_{x}^{2}\phi)$ has scaling
dimension

$\displaystyle\Delta_{\Delta_{x}^{2}\phi}$ $\displaystyle=\frac{4}{\pi
R^{2}}\int_{0}^{\pi}dx\,\sin^{3}(x/2)=\frac{16}{3\pi R^{2}},$ (36)

and as such is only irrelevant provided $R^{2}<8/3\pi$. Since $3\pi/8>1$,
there is still no choice of constant $R^{2}$ for which both
$\cos(\Delta_{x}^{2}\phi),\,\cos(\Delta_{x}^{2}\theta)$ are irrelevant.

Figure 5:  The same as in fig. 2, but excluding $\cos(\Phi_{q})$ operators
which have nonzero dipole moment.

For the choice of $R^{2}(k_{x},0)$ in (35), the analogue of fig. 2 in the
presence of dipole symmetry is shown in fig. 5. We see that imposing dipole
symmetry leads to a slight increase in the size of the stability region,
mostly in the region where $\lambda_{2}$ is close to $1$.

Finally, note that the analysis for the case with dipole conservation but
without translation symmetry is the same as the case with translation symmetry
but without dipole symmetry, just with $R^{2}(k_{x},0)\leftrightarrow
R^{-2}(k_{x},0)$ (as translations act as a dipole symmetry on the $\theta$
fields).

###  IV.5 Stability in the presence of large marginal deformations

We found above that operators preserving the subsystem symmetries are always
either marginal or irrelevant. For those operators which are marginal, we can
then ask whether or not our conclusions about stability change if such
operators are explicitly included in the fixed point action.

Focusing on operators built out of $\phi$, we found above that
$\cos(\Delta_{x}^{n}\Delta_{y}^{m}\phi)$ is marginal provided that
$\min(n,m)=1$. We may therefore consider adding these terms to the action and
then expanding the cosines to quadratic order, yielding an action of the form

$\displaystyle S$ $\displaystyle=\frac{R^{2}}{{2\pi}}\int
d\tau\sum_{i}\Big{(}\frac{1}{2\varsigma}(\partial_{\tau}\phi)^{2}+\varsigma(\Delta_{x}\Delta_{y}\phi)^{2}+\varsigma\sum_{n=2}^{\infty}\left(\alpha_{x,n}(\Delta_{x}^{n}\Delta_{y}\phi)^{2}+\alpha_{y,n}(\Delta_{x}\Delta_{y}^{n}\phi)^{2}\right)\Big{)},$
(37)

where the $\alpha_{a,n}$ are dimensionless constants, whose magnitudes we
generically we expect to be smaller for larger values of $n$ (since
microscopically such terms correspond to $2^{n+1}$-body boson operators).
Nevertheless, as the added terms are marginal, from a field theory perspective
it makes sense to treat them on the same footing as the leading
$(\Delta_{x}\Delta_{y}\phi)^{2}$ ring-exchange term.

The terms in the second line of (37) modify the dispersion relation of $\phi$
to

$\displaystyle\varepsilon_{\mathbf{k}}$
$\displaystyle=4\varsigma|\sin(k_{x}a/2)\sin(k_{y}a/2)|\left(1+\sum_{n=2}^{\infty}4^{n-1}\left(\alpha_{x,n}\sin^{2n-2}(k_{x}a/2)+\alpha_{y,n}\sin^{2n-2}(k_{y}a/2)\right)\right)^{1/2}.$
(38)

Since the added operators preserve all subsystem symmetries, they do not
change the fact that the dispersion vanishes along the coordinate axes in
momentum space. We will assume that the $\alpha_{a,n}$ are generic enough such
that the dispersion continues to vanish only along the coordinate axes, that
the first derivatives $(\partial_{k_{x}}\varepsilon)(0,k_{y})$ and
$(\partial_{k_{y}}\varepsilon)(k_{x},0)$ continue to be non-vanishing for all
nonzero $k_{y}$ and $k_{x}$ respectively, and that
$(\partial_{k_{x}}\partial_{k_{y}}\varepsilon)({\mathbf{0}})\neq 0$. With
these assumptions the RG eigenvalues of operators will still be determined by
integrals over momentum shells surrounding the coordinate axes, and can
computed using the same analysis as that appearing below (22). To leading
order in $\eta=\Lambda a\ll 1$ (with the low energy modes defined as having
momentum satisfying $\varepsilon_{\mathbf{k}}\leq\varsigma\eta^{2}$, as
before), we find

$\displaystyle\Delta_{\Phi_{q}}$
$\displaystyle=\frac{a}{{2\pi}}\int_{k_{0}}^{\pi/a}\left(\frac{dk_{x}\,|q(k_{x},0)|^{2}}{R^{2}(k_{x},0)[(\partial_{k_{y}a}\varepsilon)(k_{x},0)]}+\frac{dk_{y}\,|q(0,k_{y})|^{2}}{R^{2}(0,k_{y})[(\partial_{k_{x}a}\varepsilon)(0,k_{y})]}\right),$
(39)

where
$k_{0}=\Lambda/[(\partial_{k_{x}a}\partial_{k_{y}a}\varepsilon)({\mathbf{0}})]$
(which may be sent to zero when computing the scaling dimensions of the dipole
operators). As a sanity check, note that this more general expression reduces
to (23) upon setting
$\varepsilon_{\mathbf{k}}=4|\sin(k_{x}a/2)\sin(k_{y}a/2)|$.

By tuning the $\alpha_{a,n}$ appropriately it may very well be possible to
render all of the dipole operators irrelevant, thereby producing an EBL phase
stabilized by $(b\geq 8)$-body ring-exchange terms. We leave a more detailed
investigation of this possibility to the future.

###  IV.6 Stability in the continuum limit

All of the analysis in this paper has been concerned with the thermodynamic
limit, where the lattice spacing is kept finite as the system size is sent to
infinity. While this is the limit relevant to doing condensed matter physics,
other types of limits, particularly continuum limits in which the lattice
spacing is sent to zero, can be more natural in field theory settings.

An interesting aspect of the EBL is that the physics is sensitive to the type
of limit which is taken (see e.g. Seiberg and Shao (2020a); Gorantla et al.
(2021a)). In particular, in the continuum limit taken in Ref. Gorantla et al.
(2021a), it was shown that all operators violating the particle and vortex
subsystem symmetries are infinitely irrelevant, in the sense of having
ultralocal correlation functions Gorantla et al. (2021a). Note that this
includes the dipolar operators that we found to be responsible for
destabilizing the $R^{2}={\rm const.}$ fixed point in our calculations. This
means that conclusions about stability in the EBL depend crucially on the type
of limit which is taken. It would be interesting to understand at a deeper
level why this is so.

##  V Generalization to 3d

Everything we have discussed so far admits a straightforward generalization to
bosons hopping on a cubic lattice in 3+1D, where we may consider a model whose
kinetic term is dominated by a cube-exchange term , as studied in refs You et
al. (2020b); Seiberg and Shao (2020b); Gorantla et al. (2020). The appropriate
analogue of the free action (6) is

$S=\frac{R^{2}}{{2\pi}}\int
d\tau\sum_{i}\left(\frac{1}{2\varsigma}(\partial_{\tau}\phi)^{2}+\varsigma(1-\cos(\Delta_{x}\Delta_{y}\Delta_{z}\phi))\right),$
(40)

where $\phi$ is again a field which keeps track of the phase of the UV bosons,
with the boson density being written in terms of a dual field $\theta$ as
$n=\mkern 1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mu+\frac{1}{{2\pi}a^{3}}\Delta_{x}\Delta_{y}\Delta_{z}\theta$. The action
of translation symmetry on $\theta$ is, in analogy to (4),

$T_{\bm{\mu}}:\,\theta(\mathbf{r})\mapsto\theta(\mathbf{r}+{\bm{\mu}})+{2\pi}\mkern
1.5mu\overline{\mkern-1.5mun\mkern-1.5mu}\mkern
1.5mu\left((x+\mu_{x})(y+\mu_{y})(z+\mu_{z})-xyz\right).$ (41)

Again as in the 2+1D case, this model as written has an infinite-dimensional
group of linear subsystem symmetries, with the boson number being separately
conserved along every line parallel to one of the coordinate axes. Again as in
the 2+1D case, part of our task is to determine whether or not terms which
break this gigantic symmetry group are relevant (in the technical sense).

If we as before generalize to let $R^{2}$ be a function of $\mathbf{k}$ and
work in the phase where the cube-exchange term dominates so that the cosine
may be Taylor expanded, we obtain the Gaussian action

$\displaystyle
S_{0}=\int_{\mathbf{k},\omega}\frac{R^{2}(\mathbf{k})}{{4\pi}\varsigma}\left(\omega^{2}+\varepsilon_{\mathbf{k}}^{2}\right)|\phi(\omega,\mathbf{k})|^{2},$
(42)

where we have written
$\int_{\mathbf{k},\omega}=\int\frac{d^{2}k}{({2\pi}/a)^{3}}\frac{d\omega}{{2\pi}}$,
and where the dispersion is now

$\varepsilon_{\mathbf{k}}=8\varsigma|\sin(k_{x}a/2)\sin(k_{y}a/2)\sin(k_{z}a/2)|.$
(43)

The dispersion (43) vanishes along the codimension-1 surface in momentum space
spanned by the planes where one of $k_{x},k_{y},k_{z}$ vanishes. Given the
analysis of the preceding sections it should be clear how to set up RG:
scaling dimensions are controlled by the function $R^{2}(k_{x},k_{y},0)$
(assuming cubic symmetry), and the RG proceeds by integrating out shells with
momentum satisfying

$(\eta^{\prime})^{3}<\varepsilon_{\mathbf{k}}/\varsigma<\eta^{3},$ (44)

with $\eta=\Lambda a\ll 1$ and $\eta^{\prime}=\eta(1-dt)$ as before.

We now define the scaling dimension of an operator $\cos(\mathcal{O})$ in
terms of the associated fast mode propagator as (cf (18))

$G_{>}(0,{\mathbf{0}})=6\Delta_{\mathcal{O}}dt+O(dt^{2}).$ (45)

The factor of 6 on the RHS is chosen so that correlation functions of
$\mathcal{O}$ are functions of spacetime distances to the power of
$2\Delta_{\mathcal{O}}$ (as can be shown along the lines of the calculations
in appendix A). The RG eigenvalue of a coupling $g$ associated with
$\cos(\mathcal{O})$ is consequently

$y_{g}=d_{\mathcal{O}}-3\Delta_{\mathcal{O}},$ (46)

with $d_{\mathcal{O}}$ determined as before by requiring that when evaluated
on typical field configurations, $\eta^{d_{\mathcal{O}}}$ times the
perturbation goes as $(\Delta_{x}\Delta_{y}\Delta_{z}\phi)^{2}\sim\eta^{6}$
(for example if $\mathcal{O}=\Delta_{x}\Delta_{y}\phi$, then
$d_{\mathcal{O}}=6$).

Evaluating $G_{>}$, we see that the scaling dimension of a general operator
$\cos(\Phi_{q})$ is

$\displaystyle\Delta_{\Phi_{q}}$ $\displaystyle=\frac{1}{6dt\cdot
8\pi^{2}}\int_{dS_{\Lambda}}\frac{dk_{x}\,dk_{y}\,dk_{z}}{R^{2}(\mathbf{k})\Lambda^{3}}|q(\mathbf{k})|^{2}$
(47)
$\displaystyle=\frac{1}{32\pi^{2}}\left(\int_{-\pi/a}^{-\Lambda}+\int_{\Lambda}^{\pi/a}\right)dk_{y}\,dk_{z}\,$
$\displaystyle\qquad\times\frac{|q(0,k_{y},k_{z})|^{2}}{R^{2}(0,k_{y},k_{z})|\sin(k_{y}a/2)\sin(k_{z}a/2)|}+\cdots,$

where $\cdots$ is a stand-in for analogous integrals over $dk_{x}\ dk_{z}$ and
$dk_{x}\ dk_{y}$. From (47) we again see that $\Delta_{\Phi_{q}}$ diverges as
$\eta\rightarrow 0$ unless $q({\mathbf{0}})=0$, i.e. unless $\cos(\Phi_{q})$
is neutral under the global $U(1)$ boson number conservation. However unlike
the two-dimensional case, we see that $\Delta_{\Phi_{q}}$ also diverges unless
$q(\mathbf{k})$ vanishes when any two of $k_{x},k_{y},k_{z}$ vanish. Any
operator with finite scaling dimension must therefore be neutral under all
planar subsystem symmetries, i.e. must separately conserve the number of
bosons in each lattice plane. Finally, we see that if $q(\mathbf{k})=0$
whenever any one of $k_{x,y,z}=0$—i.e. if ${\cos(\Phi_{q})}$ is neutral under
the linear subsystem symmetries—we have $\Delta_{\Phi_{q}}=0$. As in the 2+1D
case, the latter type of operators have $d_{\mathcal{O}}\leq 0$, and as such
they are always either marginal or irrelevant. All of the preceding statements
apply equally well to $\cos(\Theta_{q})$ operators, with the only change being
$R^{2}(\mathbf{k})\leftrightarrow R^{-2}(\mathbf{k})$ in (47).

From the above discussion, as long as we are only interested in operators with
the potential to destabilize the Gaussian fixed point, we may without loss of
generality restrict our attention to operators which are invariant under the
planar subsystem symmetries and which have e.g.
$q(0,k_{y},k_{z})=q(k_{x},0,k_{z})=0,q(k_{x},k_{y},0)\neq 0$ (i.e., operators
which have vanishing dipole moment, and have nonzero quadrupole moment
oriented along ${\mathbf{\hat{z}}}$). Any $\Phi_{q}$ fitting the bill may be
written as

$\displaystyle\Phi_{q}$
$\displaystyle=\sum_{i}q_{i}\big{[}\phi(\mathbf{r}_{i})-\phi(\mathbf{r}_{i}+{\mathbf{\hat{x}}}a)+\phi(\mathbf{r}_{i}+{\mathbf{\hat{x}}}a+{\mathbf{\hat{y}}}a)$
(48) $\displaystyle\qquad-\phi(\mathbf{r}_{i}+{\mathbf{\hat{y}}}a)\big{]},$

where the $\\{q_{i}\\}$ are integers and where the sum is over the sites of a
two-dimensional square lattice. In terms of the $\\{q_{i}\\}$, the scaling
dimensions of interest may then be written as

$\displaystyle\Delta_{\Phi_{q}}$
$\displaystyle=\frac{1}{2\pi^{2}}\int_{-\pi/a}^{\pi/a}dk_{x}\,dk_{y}\frac{|\sin(k_{x}a/2)\sin(k_{y}a/2)|}{R^{2}(k_{x},k_{y},0)}$
(49)
$\displaystyle\qquad\times\sum_{i,j}q_{i}q_{j}\cos([\mathbf{r}_{i}-\mathbf{r}_{j}]\cdot\mathbf{k}),$
$\displaystyle\Delta_{\Theta_{q}}$
$\displaystyle=\frac{1}{2\pi^{2}}\int_{-\pi/a}^{\pi/a}dk_{x}\,dk_{y}\,|\sin(k_{x}a/2)\sin(k_{y}a/2)|$
$\displaystyle\qquad\times
R^{2}(k_{x},k_{y},0)\sum_{i,j}q_{i}q_{j}\cos([\mathbf{r}_{i}-\mathbf{r}_{j}]\cdot\mathbf{k}).$

These operators all have $d_{\mathcal{O}}=6$, and hence from (46) their
relevance is determined by comparing $\Delta_{\mathcal{O}}$ with 2. Also note
that by (41), $\cos(\Theta_{q})$ is translation-invariant only if its (vortex)
monopole, dipole, and quadrupole moments vanish. Therefore translation
symmetry restricts to $\\{q_{i}\\}$ such that $\sum_{i}q_{i}=0$ in the second
line above.

####  V.0.1 Constant $R^{2}$

Consider first the case where $R^{2}$ is independent of $\mathbf{k}$. The
simplest potentially relevant cosine of the $\phi$ variables is
$\cos(\Delta_{x}\Delta_{y}\phi)$, which has scaling dimension

$\Delta_{\Delta_{x}\Delta_{y}\phi}=\frac{8}{\pi^{2}R^{2}},$ (50)

therefore being irrelevant only when $R^{2}<4/\pi^{2}\approx 0.4$. On the
other hand, the simplest potentially relevant translation-invariant cosine of
the $\theta$ variables is e.g. $\cos(\Delta_{x}^{2}\Delta_{y}\theta)$, with

$\Delta_{\Delta_{x}^{2}\Delta_{y}\theta}=\frac{64R^{2}}{3\pi^{2}},$ (51)

which is irrelevant only if $R^{2}>3\pi^{2}/32\approx 0.9$. Therefore the
theory with constant $R^{2}$ is always unstable, as in the 2+1D case. However,
if one imposes a global quadrupole symmetry on the $\phi$ fields, the simplest
allowed $\cos(\Phi_{q})$ operator is then
$\Phi_{q}=\Delta_{x}^{2}\Delta_{y}\phi$, which is irrelevant provided that
$R^{2}<32/(3\pi^{2})\approx 1.1$. There is then a small region $0.9\lessapprox
R^{2}\lessapprox 1.1$ for which both this operator and
$\cos(\Delta_{x}^{2}\Delta_{y}\theta)$ are irrelevant.

Figure 6:  Stability regions for the XY-cube model with the choice of
$R^{2}(\mathbf{k})$ given in (52); translation symmetry is assumed.

####  V.0.2 General $R^{2}(\mathbf{k})$

For general choices of $R^{2}(\mathbf{k})$ the story is similar to the 2d
case, expect with slightly larger regions of stability for choices of
$R^{2}(\mathbf{k})$ analogous to that of (35). An example of the stability
region in the $\lambda_{1}$-$\lambda_{2}$ plane for the choice

$R(\mathbf{k})=\lambda_{1}(1+\lambda_{2}\cos(a[k_{x}+k_{y}+k_{z}]))^{2}$ (52)

is shown in fig. 6.

##  VI Conclusion

In this note we have discussed a natural scheme for performing RG in the
exciton Bose liquid and related models. We showed that although the simplest
type of exciton Bose liquid is unstable within our RG scheme, a certain choice
of marginal deformations can be made such that a stable phase is likely to be
realizable. This last point was argued for on the basis of a simple numerical
search, and it would be nice to obtain an analytic perspective on this issue,
perhaps along the lines of that developed in Plamadeala et al. (2014).

More generally, this way of thinking about RG in models whose IR fixed points
involve the appearance of a microscopic length scale may be useful in other
contexts, e.g. in studying the 3+1D XY plaquette model Seiberg and Shao
(2020b); Gorantla et al. (2021a). It would also be interesting to study RG
flows in these models more generally, beyond just the rather elementary
evaluation of RG eigenvalues performed here.

## Acknowledgments

I thank the University of Colorado Boulder and Marvin Qi for hospitality while
this work was written up, and am grateful to Ho-Tat Lam and Shu-Heng Shao for
discussions and detailed feedback on a first draft. I am supported by the
Hertz Fellowship.

##  Appendix A Correlation functions

In this appendix we will compute a few correlation functions at the 2+1D EBL
fixed point (the appropriate generalization to the 3+1D example of section V
is straightforward). See Paramekanti et al. (2002); Gorantla et al. (2021a)
for a detailed analysis of related correlation functions in various different
limits.

We will focus on correlation functions of exponentials of $\Phi_{q}$ operators
defined as in (21); correlators involving the $\theta$ fields can be obtained
by inverting $R^{2}(\mathbf{k})$, as usual. Letting $S_{\Lambda}$ denote the
low-energy region in momentum space (where
$\varepsilon_{\mathbf{k}}\leq\varsigma\eta^{2}$), the two-point function of
$e^{i\Phi_{q}}$ is (working in units where $\varsigma=1$ for simplicity)

$\displaystyle-\ln\langle
e^{i\Phi_{q}(\mathbf{r},\tau)}e^{-i\Phi_{q}({\mathbf{0}},0)}\rangle$
$\displaystyle=\frac{a^{2}}{{4\pi}}\int_{S_{\Lambda}}\frac{dk_{x}\,dk_{y}}{R^{2}(\mathbf{k})}\frac{|q(\mathbf{k})|^{2}}{4|\sin(k_{x}a/2)\sin(k_{y}a/2)|}$
(53)
$\displaystyle\qquad\times\left(1-(\cos(xk_{x})\cos(yk_{y})-\sin(xk_{x})\sin(yk_{y}))e^{-4\tau|\sin(k_{x}a/2)\sin(k_{y}a/2)|}\right)$

where $\mathbf{r}=(x,y)$. Consider what needs to happen in order to cancel the
logarithmic divergences that could arise when $k_{x},k_{y}$ are small. Sending
$k_{x}\rightarrow 0$ tells us that the RHS diverges unless either $y=0$ or
$q(0,k_{y})=0\,\,\forall\,\,k_{y}$, and likewise with $x\leftrightarrow y$.
Therefore if $\mathbf{r}\neq 0$ is parallel to ${\mathbf{\hat{x}}}$ (to
${\mathbf{\hat{y}}}$), the correlation function vanishes in the thermodynamic
limit unless $\Phi_{q}$ individually conserves the number of bosons along each
column (each row) of the lattice. If $\mathbf{r}\neq 0$ is not parallel to
either of the coordinate axes, the correlation function vanishes unless
$\Phi_{q}$ respects both subsystem symmetries. In fact in this case, the
correlator is asymptotically constant. Indeed, at all points in $S_{\Lambda}$,
we may always write $q(\mathbf{k})$ as either $q(k_{x},0)$ or $q(0,k_{y})$, up
to corrections vanishing as $\eta\rightarrow 0$. As such, if
$q(k_{x},0)=q(0,k_{y})=0$ for all $k_{x},k_{y}$, the correlator is constant up
to terms that vanish with $\eta\rightarrow 0$.

Consider now for simplicity the equal-time correlator, with $\tau=0$. From the
comments above, the only interesting case is one where e.g. $\mathbf{r}=(0,y)$
and $q(0,k_{y})=0\,\,\forall\,\,k_{y}$, but where $q(k_{x},0)$ is nontrivial,
and $y\neq 0$. Then working up to terms that vanish as $\eta\rightarrow 0$, we
have

$\displaystyle-\ln\langle
e^{i\Phi_{q}(\mathbf{r},\tau)}e^{-i\Phi_{q}({\mathbf{0}},0)}\rangle$
$\displaystyle=\frac{a^{2}}{{2\pi}}\int_{\Lambda}^{\pi/a}dk_{x}\int_{0}^{\eta^{2}/2a\sin(k_{x}a/2)}dk_{y}\frac{|q(k_{x},0)|^{2}(1-\cos(yk_{y}))}{R^{2}(k_{x},0)\sin(k_{x}a/2)k_{y}a}.$
(54)

We now consider the limit of large spatial separation, where $y/a\gg
1/\eta^{2}$. In this regime we may perform the integral over $k_{y}$ to give

$\displaystyle-\ln\langle
e^{i\Phi_{q}(\mathbf{r},\tau)}e^{-i\Phi_{q}({\mathbf{0}},0)}\rangle$
$\displaystyle=\frac{a}{{2\pi}}\int_{\Lambda}^{\pi/a}\frac{dk_{x}\,|q(k_{x},0)|^{2}}{R^{2}(k_{x},0)\sin(k_{x}a/2)}\ln(y\eta^{2}/a),$
(55)

where we have only kept the leading piece in the $y/a\rightarrow\infty$ limit.
This means that for large $y$, we have

$\displaystyle\langle
e^{i\Phi_{q}(0,y,0)}e^{-i\Phi_{q}({\mathbf{0}},0)}\rangle\sim\frac{1}{|ya\Lambda^{2}|^{2\Delta_{\Phi_{q}}}},$
(56)

where $\Delta_{\Phi_{q}}$ is the scaling dimension as determined by taking
$q(0,k_{y})=0$ in (23). Note that as claimed in the main text, our definition
of $\Delta_{\mathcal{O}}$ ensures that the correlator is proportional to
$|y|^{-2\Delta_{\Phi_{q}}}$.

On the other hand, consider the case where $\mathbf{r}=0,\tau\neq 0$. We then
have

$\displaystyle-\ln\langle
e^{i\Phi_{q}({\mathbf{0}},\tau)}e^{-i\Phi_{q}({\mathbf{0}},0)}\rangle$
$\displaystyle=\frac{a^{2}}{{4\pi}}\int_{S_{\Lambda}}\frac{dk_{x}\,dk_{y}}{R^{2}(\mathbf{k})}\frac{|q(\mathbf{k})|^{2}(1-e^{-4\tau|\sin(k_{x}a/2)\sin(k_{y}a/2)|})}{4|\sin(k_{x}a/2)\sin(k_{y}a/2)|}$
(57)
$\displaystyle=\frac{a}{2\pi}\int_{\Lambda}^{\pi/a}dk_{x}\int_{0}^{\eta^{2}/2a\sin(k_{x}a/2)}dk_{y}\frac{|q(k_{x},0)|^{2}(1-e^{-2\tau
k_{y}a\sin(k_{x}a/2)})}{\sin(k_{x}a/2)k_{y}}+(k_{x}\leftrightarrow k_{y}).$

In the large $\tau$ limit where $\tau\gg 1/\eta^{2}$ (recall that we are in
units where $\varsigma=1$) we may do the integrals whose upper limits go as
$\eta^{2}$, and one can check that we obtain

$\displaystyle\langle
e^{i\Phi_{q}({\mathbf{0}},\tau)}e^{-i\Phi_{q}({\mathbf{0}},0)}\rangle\sim\frac{1}{|\tau\Lambda^{2}|^{2\Delta_{\Phi_{q}}}},$
(58)

where $\Delta_{\Phi_{q}}$ is again as in (23). The correlation functions when
both $\tau$ and $\mathbf{r}$ are nonzero are obtained in a similar way.

## References

  * Paramekanti et al. (2002) A. Paramekanti, L. Balents, and M. P. Fisher, Physical Review B 66, 054526 (2002).
  * You et al. (2020a) Y. You, J. Bibo, F. Pollmann, and T. L. Hughes, arXiv preprint arXiv:2008.01746 (2020a).
  * You et al. (2020b) Y. You, Z. Bi, and M. Pretko, Physical Review Research 2, 013162 (2020b).
  * Gorantla et al. (2021a) P. Gorantla, H. T. Lam, N. Seiberg, and S.-H. Shao, arXiv preprint arXiv:2108.00020 (2021a).
  * Seiberg and Shao (2020a) N. Seiberg and S. Shao, arXiv preprint arXiv:2003.10466 (2020a).
  * Tay et al. (2011) T. Tay, O. I. Motrunich, et al., Physical Review B 83, 205107 (2011).
  * Kadanoff (1966) L. P. Kadanoff, Physics Physique Fizika 2, 263 (1966).
  * You and Moessner (2021) Y. You and R. Moessner, arXiv preprint arXiv:2106.07664 (2021).
  * You et al. (2021) Y. You, J. Bibo, T. L. Hughes, and F. Pollmann, arXiv preprint arXiv:2101.01724 (2021).
  * Gorantla et al. (2021b) P. Gorantla, H. T. Lam, N. Seiberg, and S.-H. Shao, arXiv preprint arXiv:2103.01257 (2021b).
  * Shankar (1994) R. Shankar, Reviews of Modern Physics 66, 129 (1994).
  * Lake et al. (2021) E. Lake, T. Senthil, and A. Vishwanath, Physical Review B 104, 014517 (2021).
  * Lake (to appear) E. Lake (to appear).
  * Vishwanath et al. (2004) A. Vishwanath, L. Balents, and T. Senthil, Physical Review B 69, 224416 (2004).
  * Fradkin et al. (2004) E. Fradkin, D. A. Huse, R. Moessner, V. Oganesyan, and S. L. Sondhi, Physical Review B 69, 224415 (2004).
  * Lam (private communication) H.-T. Lam (private communication).
  * Xu and Fisher (2007) C. Xu and M. P. Fisher, Physical Review B 75, 104428 (2007).
  * Plamadeala et al. (2014) E. Plamadeala, M. Mulligan, and C. Nayak, Physical Review B 90, 241101 (2014).
  * Lake and Hermele (2021) E. Lake and M. Hermele, arXiv preprint arXiv:2107.09073 (2021).
  * Mukhopadhyay et al. (2001) R. Mukhopadhyay, C. Kane, and T. Lubensky, Physical Review B 64, 045120 (2001).
  * Vishwanath and Carpentier (2001) A. Vishwanath and D. Carpentier, Physical review letters 86, 676 (2001).
  * Seiberg and Shao (2020b) N. Seiberg and S.-H. Shao, arXiv preprint arXiv:2004.00015 (2020b).
  * Gorantla et al. (2020) P. Gorantla, H. T. Lam, N. Seiberg, and S.-H. Shao, SciPost Phys 9, 2007 (2020).