# Anisotropic deconfined criticality in Dirac spin liquids

Henry Shackleton  Department of Physics, Harvard University, Cambridge MA
02138, USA  Subir Sachdev  Department of Physics, Harvard University,
Cambridge MA 02138, USA  School of Natural Sciences, Institute for Advanced
Study, Princeton NJ 08540, USA

###### Abstract

We analyze a Higgs transition from a $\operatorname{U}(1)$ Dirac spin liquid
to a gapless $\mathbb{Z}_{2}$ spin liquid. This $\mathbb{Z}_{2}$ spin liquid
is of relevance to the spin $S=1/2$ square lattice antiferromagnet, where
recent numerical studies have given evidence for such a phase existing in the
regime of high frustration between nearest neighbor and next-nearest neighbor
antiferromagnetic interactions (the $J_{1}$-$J_{2}$ model), appearing in a
parameter regime between the vanishing of Néel order and the onset of valence
bond solid ordering. The proximate Dirac spin liquid is unstable to monopole
proliferation on the square lattice, ultimately leading to Néel or valence
bond solid ordering. As such, we conjecture that this Higgs transition
describes the critical theory separating the gapless $\mathbb{Z}_{2}$ spin
liquid of the $J_{1}$-$J_{2}$ model from one of the two proximate ordered
phases. The transition into the other ordered phase can be described in a
unified manner via a transition into an unstable $\operatorname{SU}(2)$ spin
liquid, which we have analyzed in prior work. By studying the deconfined
critical theory separating the $\operatorname{U}(1)$ Dirac spin liquid from
the gapless $\mathbb{Z}_{2}$ spin liquid in a $1/N_{f}$ expansion, with
$N_{f}$ proportional to the number of fermions, we find a stable fixed point
with an anisotropic spinon dispersion and a dynamical critical exponent $z\neq
1$. We analyze the consequences of this anisotropic dispersion by calculating
the angular profiles of the equal-time Néel and valence bond solid correlation
functions, and we find them to be distinct. We also note the influence of the
anisotropy on the scaling dimension of monopoles.

###### Contents

  1. I Introduction
  2. II Summary of prior work and derivation of continuum theory
  3. III Large $N_{f}$ effective action
  4. IV Renormalization group analysis
    1. IV.1 Fermion self-energy
    2. IV.2 Fixed points
    3. IV.3 Néel and VBS order parameter corrections
      1. IV.3.1 One-loop vertex corrections
    4. IV.4 Tree-level effect of velocity anisotropy on correlation functions
  5. V Monopoles
  6. VI Conclusions and future directions
  7. VII Acknowledgements
  8. A Alternative proximate $\mathbb{Z}_{2}$ spin liquids
    1. A.1 Symmetry fractionalization in the continuum staggered flux theory
    2. A.2 Translationally-invariant spin liquid ansatzes
    3. A.3 Non-translationally-invariant spin liquid ansatzes
  9. B Emergent subsystem symmetries in the $\operatorname{SU}(2)\rightarrow\mathbb{Z}_{2}$ transition
  10. C Calculation of effective bosonic propagators
  11. D Derivation of one-loop renormalization group equations
  12. E Anisotropic correlation functions in real space
  13. F Perturbative corrections to monopole scaling dimension

##  I Introduction

Quantum antiferromagnetism has been a topic of intense study for many decades,
and has led to many new insights into quantum many-body phenomenon. In
particular, a new class of quantum phases, known as quantum spin liquids
(QSLs) [1, 2, 3], are predicted to emerge in certain parameter regimes of
antiferromagnets due to the combination of geometric frustration and quantum
fluctuations. A generic feature of these QSLs is the existence of
fractionalized excitations, which cannot be created individually by local
operators.

A particularly well-studied antiferromagnetic model is the $J_{1}$-$J_{2}$
model on the square lattice [4, 5, 6, 7, 8, 9], which has nearest-neighbor and
next-nearest-neighbor antiferromagnetic exchange interactions with
coefficients $J_{1}$ and $J_{2}$, respectively. It is known that the ground
state of the nearest-neighbor square lattice antiferromagnet ($J_{2}=0$) has
long-range Néel order i.e. global $\operatorname{SU}(2)$ spin rotation
symmetry is broken with the spin expectation value
$\langle{\bm{S}}_{{\bm{i}}}\rangle=\eta_{\bm{i}}{\bm{N}}_{0}$ , where
${\bm{S}}_{\bm{i}}$ is the spin operator on site ${\bm{i}}$,
$\eta_{\bm{i}}=\pm 1$ on the two checkerboard sublattices, and ${\bm{N}}_{0}$
is the antiferromagnetic moment. The next-nearest-neighbor antiferromagnetic
interactions compete against this Néel order, and the nature of this model as
a function of $J_{2}/J_{1}$, in particular in the regime of high frustration,
$J_{2}/J_{1}\approx 0.5$, remain a key open question.

An early proposal [10, 11, 8, 9] was that there was a direct transition from
the Néel state to a valence bond solid (VBS) (see Fig. 1) which restores spin
rotation symmetry but breaks lattice symmetries. This is followed by a first
order transition at larger $J_{2}/J_{1}$ to a ‘columnar’ state which breaks
spin rotation symmetry, and which we do not address in the present work. This
led to the development of a theory of ‘deconfined criticality’ [12, 13, 14]
which allowed for a direct transition between two symmetry-broken phases, a
phenomenon disallowed by conventional Landau-Ginzberg theory. Numerical
evidence has since accumulated for the presence of a VBS phase in the
$J_{1}$-$J_{2}$ model - in particular, it was recently shown that a non-zero
antiferromagnetic third-nearest-neighbor interaction $J_{3}$ stabilizes a
clear VBS phase in a large parameter range of $J_{2}/J_{1}$ [15], and this
phase is argued to be stable down to $J_{3}=0$. The nature of the Néel-VBS
transition in this model has remained a question of significant debate.
However, in the past two years, a consensus appears to have emerged [16] among
groups investigating this question by different numerical methods [17, 18, 19,
20], and is summarized in Fig. 1: there is a narrow window with a gapless spin
liquid phase between the Néel and VBS states. This gapless phase has been
identified [21, 22, 23, 24, 18] as a $\mathbb{Z}_{2}$ spin liquid [8, 9, 25,
26] with gapless, fermionic, $S=1/2$ spinon excitations with a Dirac-like
dispersion [27, 28, 22, 29, 30]. Although there is less of a consensus over
the precise nature of this spin liquid, variational wavefunction studies using
a Gutzwiller projection [24, 18] identify the ground state as corresponding to
the spin liquid Z2A$zz13$ in Wen’s classification [28].

In a recent work with A. Thomson [31], we proposed a unified theory that
contains multiple instabilities of the gapless $\mathbb{Z}_{2}$ spin liquid
Z2A$zz13$, which we conjecture correspond to both the neighboring VBS and Néel
orders shown in Fig. 1. The mechanism for these instabilities consists of
considering the $\mathbb{Z}_{2}$ spin liquid as a condensed phase of a parent
$\operatorname{SU}(2)$ gauge theory coupled to Higgs bosons. This
$\operatorname{SU}(2)$ theory, identified as the $\pi$-flux spin liquid, also
has a proximate phase with $\operatorname{U}(1)$ gauge symmetry, known as the
staggered flux or Dirac spin liquid. Both these phases are conjectured to be
unstable on the square lattice and ultimately lead to ordered phases, which
gives us a route for explaining the Néel and VBS ordered phases predicted to
exist alongside the $\mathbb{Z}_{2}$ spin liquid in the $J_{1}$-$J_{2}$ model.
In our previous work [31], we studied the transition between the
$\operatorname{SU}(2)$ and $\mathbb{Z}_{2}$ gauge theory using field-theoretic
techniques; in this work, we complete our study by an analysis of the
$\operatorname{U}(1)$ to $\mathbb{Z}_{2}$ transition.

The contents of this work are summarized as follows. In Section II, we review
the motivation behind our continuum model and explain its connection to the
microscopic lattice theory. In Section III, we derive a large-$N_{f}$
effective theory for the $\operatorname{U}(1)$ to $\mathbb{Z}_{2}$ transition,
with $N_{f}$ the number of fermion flavors. This allows us to study the
critical theory in a systematic ${1}/{N_{f}}$ expansion. A renormalization-
group analysis is performed in Section IV, where we extract critical exponents
of the theory to leading order in ${1}/{N_{f}}$.

Figure 1: Phases of the $S=1/2$ $J_{1}$-$J_{2}$ antiferromagnet on the square
lattice, from the numerical results of Refs. [17, 18, 19, 20], all of which
agree that the spin liquid is gapless. Each ellipse in the valence bond solid
(VBS) represents a singlet pair of electrons. Lower part of figure adapted
from Ref. [16].

The structure of our critical theory bears resemblence to prior studies of
transitions between Dirac spin liquids and a gapped $\mathbb{Z}_{2}$ gauge
theory [32, 33]. The difference between this theory and ours is reflected in
different forms of Yukawa couplings between the fermions and the Higgs fields.
This difference turns out to drastically change the qualitative features of
the critical theory. The most notable difference is the lack of Lorentz
invariance in our theory. It is known [34] that the symmetries of the square
lattice permit a velocity anisotropy term in the fermion action. In the
absence of additional gapless degrees of freedom, it has been shown [35, 36,
34] that this anisotropy is irrelevant in a $1/N_{f}$ expansion. The Yukawa
couplings of previously-studied transitions preserve Lorentz invariance.
However, we will show that the choice of Yukawa coupling necessary to realize
our specific $\mathbb{Z}_{2}$ spin liquid of interest will lead to Lorentz
symmetry breaking and a dynamical critical exponent
$z=1+0.225/N_{f}+\order{1/N_{f}^{2}}$, where $N_{f}=1$ is the physical case.

##  II Summary of prior work and derivation of continuum theory

This paper is a continuation of prior work [31], which describes multiple
possible instabilities of gapless $\mathbb{Z}_{2}$ spin liquid through a
parent $\pi$-flux phase coupled to various Higgs fields. We present a brief
summary of this derivation - further details may be found in Ref. [31].

Our starting point is the fermionic spinon theory of spin liquids, which is
derived by re-epressing the spin operators in terms of spinons
$f_{{\bm{i}}\alpha}$, $\alpha=\uparrow\,,\downarrow$ at site
${\bm{i}}=(i_{x}\,,i_{y})$ of the square lattice using

$\displaystyle\mathbf{S}_{{\bm{i}}}$
$\displaystyle=\frac{1}{2}\sum_{\alpha\,,\beta}f_{{\bm{i}}\alpha}^{\dagger}\mathbf{\sigma}_{\alpha\beta}f_{{\bm{i}}\beta}$
(1)

along with the constraints

$\displaystyle f_{{\bm{i}}\alpha}^{\dagger}f_{{\bm{i}}\alpha}=1\,,\quad
f_{{\bm{i}}\alpha}f_{{\bm{i}}\beta}\epsilon_{\alpha\beta}=0\,.$ (2)

Introducing the Nambu spinor,

$\displaystyle\psi_{{\bm{i}}}$
$\displaystyle=\begin{pmatrix}f_{{\bm{i}}\uparrow}\\\
f_{{\bm{i}}\downarrow}^{\dagger}\end{pmatrix}\,,$ (3)

and Pauli matrices $\tau^{\ell}$ which act on spinor indices, we can write a
mean-field ansatz for our Hamiltonian

$\displaystyle H$
$\displaystyle=-\sum_{{\bm{i}}{\bm{j}}}\psi_{{\bm{i}}}^{\dagger}u_{{\bm{i}}{\bm{j}}}\psi_{{\bm{j}}}\,,$
(4)

where the hoppings $u_{{\bm{i}}{\bm{j}}}$ must be determined self-
consistently. The additional degrees of freedom in this representation is
reflected by an $\operatorname{SU}(2)_{g}$ gauge symmetry, under which

$\displaystyle\operatorname{SU}(2)_{g}:\psi_{{\bm{i}}}\rightarrow
U_{g,{\bm{i}}}\psi_{{\bm{i}}}\quad\,,\quad
U_{g,{\bm{i}}}\in\operatorname{SU}(2)$ (5)

and a corresponding transformation for $u_{{\bm{i}}{\bm{j}}}$. Including gauge
fluctuations are necessary to enforce the constraint in Eq. (2), which in the
Nambu spinor variables becomes
$\psi_{{\bm{i}}}^{\dagger}\tau^{\ell}\psi_{{\bm{i}}}=0$. Different spin
liquids may be described by different mean-field ansatzes
$u_{{\bm{i}}{\bm{j}}}$, which may also spontaneously break the relevant gauge
fluctuations from $\operatorname{SU}(2)$ down to $\operatorname{U}(1)$ or
$\mathbb{Z}_{2}$.

The particular spin liquid of relevance to the $J_{1}$-$J_{2}$ model can be
labeled as Z2A$zz$13 following Wen’s classification [28]. This spin liquid,
along with two relevant proximate spin liquids, can be described by the mean-
field ansatz,

$\displaystyle u_{{\bm{i}},{\bm{i}}+\hat{x}}$
$\displaystyle=\chi\,\tau^{x}-\eta\,\tau^{y}$ (6) $\displaystyle
u_{{\bm{i}},{\bm{i}}+\hat{y}}$ $\displaystyle=\chi\,\tau^{x}+\eta\,\tau^{y}$
$\displaystyle u_{{\bm{i}},{\bm{i}}+\hat{x}+\hat{y}}$
$\displaystyle=-\gamma\,\tau^{x}$ $\displaystyle
u_{{\bm{i}},{\bm{i}}-\hat{x}+\hat{y}}$ $\displaystyle=\gamma\,\tau^{x}$

The three relevant spin liquids are:

  * •

The $\pi$-flux phase with $\operatorname{SU}(2)$ gauge symmetry corresponds to
$\chi=\eta\neq 0$ and $\gamma=0$.

  * •

The staggered flux phase with $\operatorname{U}(1)$ gauge symmetry corresponds
to $\chi\,,\eta\neq 0$, $\gamma=0$, and $\chi\neq\eta$.

  * •

The $\mathbb{Z}_{2}$ spin liquid Z2A$zz$13 is obtained from the staggered flux
phase by turning on a non-zero $\gamma$.

The dispersion relation of all three phases hosts 4 Dirac cones at low energy,
and hence, all three phases are described by $N=4$ massless Dirac fermions
minimally coupled to the corresponding gauge field. The Dirac cones of the
$\pi$-flux phase have an emergent Lorentz invariance, whereas the staggered
flux and $\mathbb{Z}_{2}$ phase have anisotropic dispersion relations on a
mean-field level (note that prior studies [35, 36, 34] show an emergent
Lorentz invariance of the staggered flux phase upon including gauge
fluctuations).

The primary mechanism for our theory of the $J_{1}$-$J_{2}$ model is the
assumption that both the $\pi$-flux and staggered flux phases on the square
lattice are ultimately unstable to ordered phases, either Néel or VBS. For the
staggered flux phase, the instability arises due to the presence of monopoles,
which are allowed by the compactness of the $\operatorname{U}(1)$ gauge
theory. The scaling dimension of these monopoles have been calculated to
second-order in a $1/N$ expansion [37, 38], with $N$ the number of fermions,
and are relevant for $N=4$. Moreover, there exists a “trivial” monopole
operator that respects the microscopic symmetries of the square lattice, and
hence is allowed by symmetry in an effective Lagrangian [39]. Proliferation of
these monopoles is conjectured to lead to ordered phases, including Néel and
VBS order [40, 41]. The $\pi$-flux phase, on the other hand, has been
conjectured to be a dual description of the DQCP separating Néel and VBS order
[42], and hence is generically unstable to these phases. A unified framework
for describing the gapless spin liquid Z2A$zz$13 as well as these two
instabilities can hence be obtained from $\text{QCD}_{3}$ with $N=2$ fermion
doublets and an $\operatorname{SU}(2)$ gauge group, and coupling it to Higgs
fields whose condensation breaks the gauge group to either
$\operatorname{U}(1)$ or $\mathbb{Z}_{2}$. The precise manner in which these
Higgs fields must couple to the Dirac fermions in order to realize the
specific spin liquids of interest can be determined by taking the continuum
limit of Eq. (6) and demanding that the condensation of the Higgs fields
modifies the Dirac dispersion relation in a manner consistent with the lattice
theory. This is the procedure carried out in Ref. [31]. These couplings may
also be determined by matching the fractionalization of the square lattice
symmetries on the microscopic level with the symmetries of the continuum
theory, which is done in Appendix A of Ref. [31]. We refer the reader to Ref.
[31] for further details of this calculation; we will simply present the
resulting Lagrangian in this work.

This procedure ultimately yields a Lagrangian consisting of four Dirac
fermions $\psi$, an $\operatorname{SU}(2)$ gauge field $A_{\mu}^{a}$, and
three three-component adjoint Higgs fields $\Phi_{1,2,3}^{a}$, $a=x,y,z$. The
$\operatorname{SU}(2)$ gauge ($\sigma$) and valley ($\mu$) Pauli matrices both
rotate between the four fermion flavors.

$\displaystyle\mathcal{L}$
$\displaystyle=\mathcal{L}_{\psi}+\mathcal{L}_{\Phi}+\mathcal{L}_{\Phi\psi}$
(7) $\displaystyle\mathcal{L}_{\psi}$
$\displaystyle=i\bar{\psi}\gamma^{\mu}\left(\partial_{\mu}-iA^{a}_{\mu}\sigma^{a}\right)\psi\,.$
$\displaystyle\mathcal{L}_{\Phi}$
$\displaystyle=\sum_{i=1}^{3}D_{\mu}\Phi_{i}^{a}D^{\mu}\Phi_{i}^{a}+V(\Phi)$
$\displaystyle\mathcal{L}_{\Phi\psi}$
$\displaystyle=\Phi_{1}^{a}\,\bar{\psi}\mu^{z}\gamma^{x}\sigma^{a}\psi+\Phi_{2}^{a}\,\bar{\psi}\mu^{x}\gamma^{y}\sigma^{a}\psi+\Phi_{3}^{a}\bar{\psi}\mu^{y}\sigma^{a}(\gamma^{x}D_{y}+\gamma^{y}D_{x})\psi$
$\displaystyle V(\Phi)$
$\displaystyle=s\left(\Phi_{1}^{a}\Phi_{1}^{a}+\Phi_{2}^{a}\Phi_{2}^{a}\right)+\widetilde{s}\,\Phi_{3}^{a}\Phi_{3}^{a}+w\,\epsilon_{abc}\,\Phi_{1}^{a}\Phi_{2}^{b}\Phi_{3}^{c}$
$\displaystyle+u\left(\Phi_{1}^{a}\Phi_{1}^{a}+\Phi_{2}^{a}\Phi_{2}^{a}\right)^{2}+\widetilde{u}\left(\Phi_{3}^{a}\Phi_{3}^{a}\right)^{2}+v_{1}\left(\Phi_{1}^{a}\Phi_{2}^{a}\right)^{2}+v_{2}\left(\Phi_{1}^{a}\Phi_{1}^{a}\right)\left(\Phi_{2}^{b}\Phi_{2}^{b}\right)$
$\displaystyle+v_{3}\left[\left(\Phi_{1}^{a}\Phi_{3}^{a}\right)^{2}+\left(\Phi_{2}^{a}\Phi_{3}^{a}\right)^{2}\right]+v_{4}\left(\Phi_{1}^{a}\Phi_{1}^{a}+\Phi_{2}^{a}\Phi_{2}^{a}\right)\left(\Phi_{3}^{b}\Phi_{3}^{b}\right)+\ldots$

We have absorbed the coefficients of the Yukawa couplings into the Higgs
fields. The general form of the Higgs potential $V(\Phi)$ is constrained by
the microscopic symmetries of the square lattice, and we present only a subset
of the possible terms. The manner in which these microscopic symmetries are
embedded in the continuum theory may be derived by starting from the original
lattice model, and these transformation properties are given in Table 1. The
action of the $\operatorname{SU}(2)$ spin rotation symmetry requires a more
careful analysis and is described below. All three Higgs fields transform
trivially under this $\operatorname{SU}(2)$ symmetry. Our choice of
representing the fermionic degrees of freedom in terms of Dirac fermions
obfuscates the full $\operatorname{SU}(2)$ spin rotation symmetry, as
rotations around the $x$ or $y$ axis involve charge conjugation. However, the
$\operatorname{U}(1)$ subgroup corresponding to rotations around the $z$ axis
is simply given by a phase shift in $\psi$, $\psi\rightarrow e^{i\theta}\psi$.
The full $\operatorname{SU}(2)$ rotation symmetry may be made explicit by
writing the theory in terms of Majorana fermions [42], but this will not be
necessary for our purposes.

| $T_{x}$ | $T_{y}$ | $P_{x}$ | $P_{y}$ | $\mathcal{T}$ | $R_{\pi/2}$
---|---|---|---|---|---|---
$\Phi_{1}^{a}$ | $-$ | $+$ | $-$ | $-$ | $-$ | $-\Phi_{2}^{a}$
$\Phi_{2}^{a}$ | $+$ | $-$ | $-$ | $-$ | $-$ | $-\Phi_{1}^{a}$
$\Phi_{3}^{a}$ | $-$ | $-$ | $+$ | $+$ | $+$ | $-$
$\psi$ | $\mu^{x}\psi$ | $\mu^{z}\psi$ | $\gamma^{x}\mu^{z}\psi$ | $-\gamma^{y}\mu^{x}\psi$ | $\gamma^{0}\mu^{y}\psi$ | $e^{i\pi\gamma^{0}/4}e^{-i\pi\mu^{y}/4}\psi$
Table 1: Listed are the microscopic symmetries of the square lattice and their
action on the continuum fields. $T_{i}$ ($P_{i}$) indicate translations
(reflections) along the $i$’th axis, $\mathcal{T}$ is time-reversal symmetry,
and $R_{\pi/2}$ is a $\pi/2$ rotation. We omit the action of
$\operatorname{SU}(2)$ spin rotation symmetry as its action on the fermions
$\psi$ is non-trivial and described in the main text. Figure 2: Mean field
phase diagram of our low energy theory obtained by minimization of the Higgs
potential in Eq. (7). Dashed (solid red) lines indicate second (first) order
transitions in mean field theory. We assume the SU(2) $\pi$-flux gauge theory
confines to a Néel state, the U(1) staggered flux gauge theory confines to a
VBS state, except at their deconfined critical boundaries to Wen’s stable,
gapless $\mathbb{Z}_{2}$ spin liquid Z2A$zz13$. The dotted blue line indicates
a possible trajectory of the square lattice antiferromagnet with increasing
$J_{2}/J_{1}$.

On a mean-field level, this theory admits three phases, illustrated in Fig. 2.
When all three Higgs fields are uncondensed, we recover the $\pi$-flux phase
with $\operatorname{SU}(2)$ gauge group. When $\Phi_{3}$ condenses, we obtain
the $\operatorname{U}(1)$ staggered flux phase. The condensation of $\Phi_{1}$
and $\Phi_{2}$ yields the gapless spin liquid Z2A$zz$13\. The masses of
$\Phi_{1}$ and $\Phi_{2}$ are fixed to be equal by the microscopic square
lattice symmetry, so both condense simultaneously. Furthermore, the symmetry-
allowed cubic term $\epsilon_{abc}\Phi_{1}^{a}\Phi_{2}^{b}\Phi_{3}^{c}$ forces
$\Phi_{3}$ to condense along with $\Phi_{1}$, $\Phi_{2}$. Our conjectured
trajectory of the $J_{1}$-$J_{2}$ model as a function of $\frac{J_{2}}{J_{1}}$
is shown by the dotted blue line in Fig. 2, where transitions into either the
staggered flux or $\pi$-flux phases drive the Néel or VBS ordering. Our theory
may also be compatible with the inclusion of an antiferromagnetic third-
nearest-neighbor $J_{3}$ term, which has been shown numerically [15] to
compete against the spin liquid phase, eventually leading to a direct Néel/VBS
transition. This phenomenon can be described in our theory by a deformation of
the ${J_{2}}/{J_{1}}$ path to a trajectory in the $\pi$-flux phase; this
deformation introduces first-order transitions near the tricritical point.

The transition from the $\pi$-flux phase to the $\mathbb{Z}_{2}$ phase was
studied in [31], and the focus of our work will be the $\operatorname{U}(1)$
to $\mathbb{Z}_{2}$ transition. Before proceeding with our analysis, we
briefly summarize the results of our study of the $\operatorname{SU}(2)$ to
$\mathbb{Z}_{2}$ transition. The primary order parameters for this theory are
the masses of $\Phi_{1,2}$, so we neglect fluctuations of $\Phi_{3}$. This
critical theory is studied in a $1/N_{f}$ expansion, with $4N_{f}$ the total
number of fermions. Due to the anisotropic couplings of the Higgs fields, the
leading-order effective propagators for $\Phi_{1,2}$ are divergent along a
one-dimensional subspace in momentum space. These lines of zero modes may be
thought of as a consequence of an emergent subsystem symmetry. More details on
this perspective are presented in Appendix B. Although these divergences are
lifted by higher-order corrections, the leading-order effective action is non-
local and more relevant at long distances than the bare Higgs kinetic term.
Performing a standard momentum-shell renormalization group study of this
leading-order theory is ill-defined and necessitates the inclusion of the
irrelevant bare Higgs kinetic terms, which become “dangerously irrelevant” due
to the singular nature of the leading-order theory. At one-loop order, these
features lead to $\log^{2}$ divergences, rather than the more standard single-
logarithm divergences found in conventional field theories. Re-exponentiating
these corrections, we predict that the universal properties of correlation
functions at criticality, rather than being power law, are instead
$r^{-\alpha}\exp\left({\beta\ln^{2}r}\right)$, with $\beta=-{6}/{\pi^{2}}$ for
the VBS correlator and $-{12}/{\pi^{2}}$ for the Néel correlator, and $\alpha$
being some non-universal coefficient.

For the remainder of our work, we study the $\operatorname{U}(1)$ to
$\mathbb{Z}_{2}$ transition, which is driven by the condensation of a charge-2
complex scalar Higgs field. Upon calculating the large-$N_{f}$ effective
action, we perform a renormalization group (RG) analysis to determine the
fixed point of our critical theory. The primary observables of interest that
we will study are the dynamical critical exponent $z$, which determines the
difference between spatial and temporal scaling, and correlations of the Néel
and VBS order parameters, which are given by fermion bilinears in our
continuum theory.

We find that the presence of the massless Higgs fields strongly modify the
critical behavior from the theory of massless QED which describes the pure
staggered flux phase without Higgs fields. The staggered flux phase normally
possesses an $\operatorname{SU}(4)$ symmetry [34], which relates various
physical order parameters including Néel and VBS. Moreover, the staggered flux
phase has an emergent Lorentz invariance at low energy, as the only velocity
anisotropy term allowed by the square lattice symmetries is irrelevant in a
$1/N_{f}$ expansion. The critical Higgs fields explicitly break the
$\operatorname{SU}(4)$ symmetry and generate a non-zero velocity anisotropy,
which further breaks the $\operatorname{U}(1)$ spatial rotation symmetry down
to the $C_{4}$ symmetry present in the microscopic model. While correlation
functions still have power law decay as in more traditional critical theories,
the angular profiles of the correlation functions are modified by the velocity
anisotropy; we calculate the tree-level effects of these modifications for the
Néel and VBS correlation functions.

Recall that in the absence of critical Higgs fields, we postulate that the
staggered flux phase is unstable to monopole proliferation. A key assumption
in our analysis of the critical theory is that monopoles are rendered
irrelevant due to the presence of critical Higgs fields and the theory may be
studied using a non-compact $\operatorname{U}(1)$ gauge field. Indeed, such an
assumption is similar to the analysis of earlier studies of deconfined
criticality between Néel and VBS orders [13, 14]. This assumption is
elaborated further in Section V, where we note that in addition to the
critical Higgs fields, the presence of a non-zero anisotropy in the fermion
dispersion relation can also affect the relevance of monopoles at criticality.
In a large-$N_{f}$ expansion, this has an $\order{N_{f}}$ effect on the
scaling dimension of the monopole, in contrast to the Higgs screening, which
is $\order{1}$.

##  III Large $N_{f}$ effective action

Our goal is to study the $\operatorname{U}(1)$ to $\mathbb{Z}_{2}$ transition
of Eq. (7). Both phases have the Higgs field $\langle\Phi_{3}^{a}\rangle\neq
0$, so we fix $\Phi_{3}^{a}=\delta_{az}\Phi$, $\Phi\neq 0$. In this theory,
the $\operatorname{SU}(2)$ gauge symmetry is broken to $\operatorname{U}(1)$,
so we only consider $A_{\mu}\equiv A_{\mu}^{z}$. It is important to include
the consequences of the cubic Higgs term, whose sign determines the form of
the low energy complex Higgs. We choose a gauge where $w\Phi<0$, and the low
energy behavior can be described in terms of a single complex Higgs

$\displaystyle\mathcal{H}$
$\displaystyle=\frac{1}{2}\left(\Phi_{1}^{x}+\Phi_{2}^{y}+i(\Phi_{1}^{y}-\Phi_{2}^{x})\right)\,.$
(8)

This field transforms as a charge 2 Higgs field under the unbroken
$\operatorname{U}(1)$ symmetry, as desired. Other linear combinations of
$\Phi_{1,2}$ are massive and can be neglected for the critical theory. As we
will clarify later, we must also assume $w<0$ in order for our transition to
be continuous - if $w$ is positive, the aforementioned massive Higgs fields
will turn out to have a negative mass at the $\order{1/N_{f}}$ fixed point,
which will lead to a first-order transition.

With these definitions, the Lagrangian that describes the
$\operatorname{U}(1)\rightarrow\mathbb{Z}_{2}$ transition is

$\displaystyle\mathcal{L}_{\text{sf}}$
$\displaystyle=\mathcal{L}_{\psi}+\mathcal{L}_{\mathcal{H}}+\mathcal{L}_{\mathcal{H}\psi}$
(9) $\displaystyle\mathcal{L}_{\psi}$
$\displaystyle=i\bar{\psi}\gamma^{\mu}D_{\mu}\psi+\Phi\,\bar{\psi}\mu^{y}\sigma^{z}\left(\gamma^{y}D_{x}+\gamma^{x}D_{y}\right)\psi\,.$
$\displaystyle\mathcal{L}_{\mathcal{H}}$
$\displaystyle=s|\mathcal{H}|^{2}+\partial_{\mu}\mathcal{H}^{*}\partial^{\mu}\mathcal{H}+u|\mathcal{H}|^{4}$
$\displaystyle\mathcal{L}_{\mathcal{H}\psi}$
$\displaystyle=\mathcal{H}\bar{\psi}\left(\mu^{z}\gamma^{x}+i\mu^{x}\gamma^{y}\right)\sigma^{-}\psi+\mathcal{H}^{\ast}\bar{\psi}\left(\mu^{z}\gamma^{x}-i\mu^{x}\gamma^{y}\right)\sigma^{+}\psi\,.$

with $D_{\mu}=\partial_{\mu}-iA_{\mu}\sigma^{z}$ and
$\sigma^{\pm}\equiv\sigma^{x}\pm i\sigma^{y}$. This Lagrangian may contain
higher-order terms, but we omit these as they will turn out to be irrelevant
in a $1/N_{f}$ expansion. The term proportional to $\Phi$ is a modification to
the Lorentz-invariant fermion propagator $1/\not{p}$ allowed by the projective
symmetry group of the staggered flux phase. In the absence of the critical
Higgs field, this velocity anisotropy has a stable fixed point value of
$\Phi=0$ [35, 36, 34]. In our case, terms of this form are spontaneously
generated at one-loop order by the critical Higgs field, and $\Phi$ acquires a
non-zero value at the fixed point.

In order to study our critical theory, we proceed in a $1/N_{f}$ expansion,
with $N_{f}$ the fermion number. Since our theory only makes sense when the
number of fermions $N$ is a multiple of $4$, we define $4N_{f}=N$; in other
words, we take $N_{f}=1$ to correspond to our physical theory. At leading
order in ${1}/{N_{f}}$, our effective bosonic action takes the form

$\displaystyle\frac{S_{b}}{N_{f}}$
$\displaystyle=\int_{k}\left[s+\Gamma(k)\right]\mathcal{H}^{*}(-k)\mathcal{H}(k)+\frac{1}{2}\Pi_{\mu\nu}(k)A_{\mu}(-k)A_{\nu}(k)$
(10)

Figure 3: The effective action for the Higgs boson (left) and
$\operatorname{U}(1)$ gauge field (right) are generated by the fermions at
leading order in a $1/N_{f}$ expansion.

where the inverse propagators $\Gamma$, $\Pi$ are generated by the one-loop
fermion diagrams shown in Fig. 3. Note that we have taken the bare Higgs mass
to scale with $N_{f}$, although we will be interested in the critical theory
where we tune the Higgs mass to zero.

To calculate the effective propagators, we need the fermion propagator, which
receives corrections to its Lorentz-invariant value of $1/\not{p}$ due to a
non-zero $\Phi$. This may be treated perturbatively in $\Phi$, but the
existence of a stable fixed point turns out to not be viewable at leading
order, so we instead proceed with a non-perturbative treatment of $\Phi$. We
include further details of this calculation in Appendix C, and cite the
results in the main text. Defining the variables

$\displaystyle k_{x,\pm}\equiv k_{x}\pm\Phi k_{y}\,,\quad k_{y,\pm}\equiv
k_{y}\pm\Phi
k_{x}\,,\quad\absolutevalue{k_{\pm}}\equiv\sqrt{k_{0}^{2}+k_{x,\pm}^{2}+k_{y,\pm}^{2}}$
(11)

the effective inverse Higgs propagator (obtained from the $\Phi$-dependent
fermion propagator) is

$\displaystyle\Gamma(k)$
$\displaystyle=\frac{1}{16N_{f}(1-\Phi^{2})}\left[\frac{k_{+}^{2}+k_{0}^{2}+2k_{x,+}k_{y,+}}{\absolutevalue{k_{+}}}+\frac{k_{-}^{2}+k_{0}^{2}-2k_{x,-}k_{y,-}}{\absolutevalue{k_{-}}}\right]\,.$
(12)

Likewise, we need the general form of the effective gauge boson propagator.
The presence of a non-zero $\Phi$ modifies the gauge coupling, and hence non-
Lorentz-invariant corrections arise both from $\Phi$-dependent modifications
to the fermion propagator as well as $\order{\Phi}$ vertices. We separate this
calculation into three pieces. The first correction comes from using the
$\order{\Phi^{0}}$ vertices, but with the full fermion propagator. This one-
loop term contributes

$\displaystyle\Pi_{\mu\nu}^{(1)}(k)$
$\displaystyle=\sum_{a=\pm}\frac{k_{a}^{2}}{8N_{f}(1-\Phi^{2})\absolutevalue{k_{a}}}\left(\delta_{\mu\nu}-\frac{k_{\mu
a}k_{\nu a}}{k_{a}^{2}}\right)\,.$ (13)

The second correction comes from using one $\order{\Phi}$ vertex, which gives
the contribution

$\displaystyle\Pi_{\mu x}^{(2)}(k)$
$\displaystyle=\Pi^{(2)}_{x\mu}(k)=\sum_{a=\pm}\frac{a\Phi
k_{a}^{2}}{8N_{f}(1-\Phi^{2})\absolutevalue{k_{a}}}\left(\delta_{\mu
y}-\frac{k_{\mu a}k_{ya}}{k_{a}^{2}}\right)$ (14) $\displaystyle\Pi_{\mu
y}^{(2)}(k)$ $\displaystyle=\Pi_{y\mu}^{(2)}(k)=\sum_{a=\pm}\frac{a\Phi
k_{a}^{2}}{8N_{f}(1-\Phi^{2})\absolutevalue{k_{a}}}\left(\delta_{\mu
x}-\frac{k_{\mu a}k_{xa}}{k_{a}^{2}}\right)\,;$

There is also an extra factor of $2$ in $\Pi^{(2)}_{xx,yy}$ due to the two
possible vertex orderings. Finally, the third correction comes from using two
$\order{\Phi}$ vertices,

$\displaystyle\Pi_{xx}^{(3)}(k)$
$\displaystyle=\sum_{a=\pm}\frac{\Phi^{2}k_{a}^{2}}{8N_{f}(1-\Phi^{2})\absolutevalue{k_{a}}}\left(1-\frac{k_{ya}^{2}}{k_{a}^{2}}\right)$
(15) $\displaystyle\Pi_{yy}^{(3)}(k)$
$\displaystyle=\sum_{a=\pm}\frac{\Phi^{2}k_{a}^{2}}{8N_{f}(1-\Phi^{2})\absolutevalue{k_{a}}}\left(1-\frac{k_{xa}^{2}}{k_{a}^{2}}\right)$
$\displaystyle\Pi_{yx}^{(3)}(k)$
$\displaystyle=\Pi^{(3)}_{xy}(k)=\sum_{a=\pm}-\frac{\Phi^{2}k_{xa}k_{ya}}{8N_{f}(1-\Phi^{2})\absolutevalue{k_{a}}}$

We verify that the combined inverse propagator
$\Pi_{\mu\nu}(k)=\sum_{i=1,2,3}\Pi^{(i)}_{\mu\nu}(k)$ annihilates the vector
$(k_{0},k_{x},k_{y})$, as required by gauge invariance. Note that
$\Pi_{\mu\nu}$ requires a gauge fixing term in order to be invertable.
Followin Ref. [34], we add the following non-local gauge fixing term to the
Lagrangian

$\displaystyle\frac{1}{4\xi\absolutevalue{k}}A_{\mu}k^{\mu}k^{\nu}A_{\nu}\,.$
(16)

All gauge-invariant observables have been checked to ensure they are
independent of the choice of $\xi$.

##  IV Renormalization group analysis

We perform a renormalization group (RG) analysis of the $\order{1/N_{f}}$
effective theory. We are interested in studying the behavior of this theory
under the rescaling

$\displaystyle k$ $\displaystyle=k^{\prime}e^{-\ell}$ (17)
$\displaystyle\omega$ $\displaystyle=\omega^{\prime}e^{-z\ell}$

We also define a rescaling of the fermion fields

$\displaystyle\psi(k,\omega)$
$\displaystyle=\psi^{\prime}(k^{\prime},\omega^{\prime})e^{\frac{\ell}{2}(2+2z-\eta_{f})}$
(18)

The Higgs and gauge fields must also be suitably rescaled, although the
anomalous dimensions of these fields will not be needed to calculate our
observables of interest. In the absence of a standard boson kinetic term at
leading order in ${1}/{N_{f}}$, we define the scaling of the boson field by
performing our RG such that the Yukawa coupling remains fixed under RG.

###  IV.1 Fermion self-energy

We first evaluate the $\order{1/N_{f}}$ contributions to the fermion self-
energy, which come from both gauge and Higgs one-loop diagrams. The self-
energy is UV divergent and requires a UV cutoff $\Lambda$. The logarithmic
derivative of the fermion self-energy with respect to this cutoff takes the
general form

$\displaystyle\Lambda\derivative{\Lambda}\Sigma(k)=C_{0}k_{0}\gamma^{0}+C_{1}(k_{x}\gamma^{x}+k_{y}\gamma^{y})+C_{2}\Phi\mu^{y}\sigma^{z}(k_{x}\gamma^{y}+k_{y}\gamma^{x})$
(19)

for constants $C_{0,1,2}$. One must verify that only these terms are generated
at one-loop order, which we have done.

In order to calculate the constants $C_{i}$, we will use the momentum-shell RG
approach outlined in Ref. [43]. The regularized one-loop contribution to the
self-energy schematically takes the form

$\displaystyle\Sigma(k)=\int\frac{\differential[3]{p}}{(2\pi)^{3}}F(p+k)G(p)\mathcal{K}\left(\frac{p^{2}}{\Lambda^{2}}\right)\mathcal{K}\left(\frac{(k+p)^{2}}{\Lambda^{2}}\right)$
(20)

where $F$ and $G$ are homogeneous functions of the three-momenta, and
$\mathcal{K}(y)$ serves as a UV cutoff with the property that
$\mathcal{K}(0)=1$ and $\mathcal{K}(y)$ falls off rapidly for large $y$. In
our calculations, we take $F$ to be the fermion propagator, and $G$ to be the
boson propagator (either Higgs or gauge), along with vertex coefficients. The
fact that $F$ and $G$ are homogeneous functions allows us to remove the
explicit dependence on $\mathcal{K}$ upon taking the logarithmic derivative
and integrating by parts. We refer to Appendix D for an explicit derivation of
this, and state the result here - the logarithmic derivative of the self-
energy takes the form

$\displaystyle\Lambda\derivative{\Lambda}\Sigma(k)$
$\displaystyle=\frac{k_{\lambda}}{8\pi^{3}N_{f}}\int_{0}^{2\pi}\differential{\phi}\int_{0}^{\pi}\sin\theta\differential{\theta}\partialderivative{F(\hat{p})}{p_{\lambda}}G(\hat{p})\,.$
(21)

where $\hat{p}\equiv(\cos\theta,\sin\theta\sin\phi,\sin\theta\cos\phi)$. The
resulting integrals in Eq. (21) are fully convergent and may be evaluated
numerically, from which we can extract the coefficients $C_{0,1,2}$.

###  IV.2 Fixed points

The RG equations for the velocity anisotropy $\Phi$ are

$\displaystyle\derivative{\Phi}{\ell}=(C_{1}-C_{2})\Phi\,.$ (22)

In the absence of the Higgs field, $\Phi$ has a stable fixed point at
$\Phi=0$. The gauge field contribution to this equation has been calculated to
leading order in $\Phi$ [34], and we verify agreement with this result.

Figure 4: An evaluation of the RG flow of $\Phi$, showing a stable fixed point
at $\Phi_{c}\approx 0.46$.

The evaluation of Eq. (22) is plotted in Fig. 4. A stable fixed point is found
at $\Phi_{c}\approx 0.45765$. At this point, the dynamical critical exponent
$z$ is given by

$\displaystyle z=1-C_{0}+C_{1}=1+\frac{0.225}{N_{f}}+\order{1/N_{f}^{2}}\,.$
(23)

Recall that when we derived this critical
$\operatorname{U}(1)\rightarrow\mathbb{Z}_{2}$ theory as a component of a
parent $\operatorname{SU}(2)$ theory, we made a gauge choice such that
$w\Phi<0$, where $w$ is the coefficient of the symmetry-allowed cubic term,
$w\epsilon_{abc}\Phi_{1}^{a}\Phi_{2}^{b}\Phi_{3}^{c}$. When
$\langle\Phi_{3}^{a}\rangle=\Phi\delta_{az}$, we can diagonalize this term to
yield a mass, $w\Phi(\mathcal{H}^{*}\mathcal{H}-\mathcal{M}^{*}\mathcal{M})$,
where $\mathcal{H}$ is the combination of $\Phi_{1,2}^{x,y}$ given in Eq.
LABEL:eq:lowEnergyHiggs and $\mathcal{M}$ is a charge-2 Higgs field of a
similar form but with $x\leftrightarrow y$. If we assume $w>0$, $\Phi<0$, then
$\mathcal{H}$ will become massless first, but the fixed-point value of $\Phi$
gives a negative mass to $\mathcal{M}$, leading to a first-order transition
driven by the condensation of $\mathcal{M}$. As a consequence, we must fix our
parent $\operatorname{SU}(2)$ theory to have $w<0$ in order to yield a
continuous transition. If we had made a gauge choice such that $w\Phi>0$, then
our theory would have been driven by the condensation of $\mathcal{M}$ rather
than $\mathcal{H}$; this still leads to the gapless $\mathbb{Z}_{2}$ spin
liquid Z2Azz13, and all gauge-invariant observables at the critical point
remain the same, although the sign of $\Phi_{c}$ changes.

###  IV.3 Néel and VBS order parameter corrections

We now calculate the vertex corrections to the Néel and VBS order parameters.
These order parameters are given by fermion bilinears and can be identified
based on the action of the microscopic square lattice symmetries on the
fermions. The VBS order parameter is given by the bilinears
$\bar{\psi}\mu^{z,x}\psi$. As mentioned previously, our particular
representation obfuscates the full $\operatorname{SU}(2)$ action of spin
rotation symmetry; however, the $\operatorname{U}(1)$ subgroup generated by
rotations around the $z$-axis is given by the global $\operatorname{U}(1)$
symmetry $\psi\rightarrow e^{i\theta}\psi$ (recall that this is not the
$\operatorname{U}(1)$ gauge symmetry, which acts as $e^{i\theta\sigma^{z}}$).
As a consequence, we focus on the $z$-component of the Néel order parameter,
which is given by $\bar{\psi}\mu^{y}\psi$. The two-point correlation functions
of these bilinears are obtained by coupling them to external sources
$J_{\text{VBS, N\'{e}el}}$ and calculating $\order{N_{f}^{-1}}$ vertex
corrections, illustrated in Fig. 5. Of note are $\order{N_{f}^{-1}}$ two-loop
corrections, the form of which were first found in Ref. [34]. These two-loop
corrections to the Néel and VBS order parameters vanish in the pure staggered
flux phase - this follows immediately from taking the trace over the internal
fermion loop and noting that $\Tr\mu^{i}=0$. We verify that these diagrams
remain zero upon the inclusion of both Higgs fields and velocity anisotropy,
although this identity is less readily apparent.

Figure 5: The $\mathcal{O}(N_{f}^{-1})$ vertex corrections which contribute to
the renormalization of the Néel and VBS order parameters. The order parameter
receives corrections at one-loop order from the Higgs fields (left) and the
gauge boson (center). An additional two-loop $\mathcal{O}(N_{f}^{-1})$
contribution (right) is possible - the diagram shown involves two intermediate
Higgs propagators, but additional diagrams with gauge propagators or mixed
gauge/Higgs propagators are possible. These diagrams vanish exactly upon
performing the trace over the internal fermion indices.

####  IV.3.1 One-loop vertex corrections

We first outline the procedure from Ref. [43] for calculating the logarithmic
corrections to the vertex functions. At zero external momenta, our one-loop
vertex corrections schematically take the form

$\displaystyle\Xi_{i}=\frac{1}{N_{f}}\int\frac{\differential[3]{p}}{(2\pi)^{3}}H_{i}(p)\mathcal{K}^{3}\left(\frac{p^{2}}{\Lambda^{2}}\right)$
(24)

where $H_{i}(p)$, $i=x,y,z$, is a homogeneous function of $p$ and illustrated
in Fig. 5. The index $i$ indicates whether the $J$ vertex includes a factor of
$\mu^{x,z}$ (VBS order parameter) or $\mu^{y}$ (Néel order parameter). Once
again, we can take the logarithmic derivative and remove the explicit cutoff
dependence, leading to the equation

$\displaystyle\Lambda\derivative{\Lambda}\Xi_{i}=\frac{1}{8\pi^{3}N_{f}}\int_{0}^{2\pi}\differential{\phi}\int_{0}^{\pi}\sin\theta\differential{\theta}H_{i}(\hat{p})\equiv
B_{i}\mu^{i}\,.$ (25)

The $B_{i}$’s are not gauge-invariant by themselves, and must be combined with
the self-energy to get a gauge-invariant quantity, which at the fixed-point
value $\Phi_{c}$ gives

$\displaystyle\eta_{\text{VBS}}$ $\displaystyle=B_{x,z}+C_{0}\approx
0.06468N_{f}^{-1}+\order{N_{f}^{-2}}\,,$ (26)
$\displaystyle\eta_{\text{N\'{e}el}}$
$\displaystyle=B_{y}+C_{0}\approx-0.01634N_{f}^{-1}+\order{N_{f}^{-2}}\,.$

The Néel and VBS correlators in momentum space have the scaling form

$\displaystyle G_{\text{N\'{e}el}}(k,\omega)$
$\displaystyle=G_{\text{N\'{e}el}}(ak,a^{z}\omega)a^{2\eta_{\text{N\'{e}el}}-1}\,,$
(27) $\displaystyle G_{\text{VBS}}(k,\omega)$
$\displaystyle=G_{\text{VBS}}(ak,a^{z}\omega)a^{2\eta_{\text{VBS}}-1}\,.$

Making a Fourier transform to real space, the equal-time Néel and VBS
correlators have the power law decay

$\displaystyle G_{\text{N\'{e}el}}(r)$
$\displaystyle\sim\frac{1}{r^{3+z-2\eta_{\text{N\'{e}el}}}}$ (28)
$\displaystyle G_{\text{VBS}}(r)$
$\displaystyle\sim\frac{1}{r^{3+z-2\eta_{\text{VBS}}}}\ $

Note that both the anomalous dimensions for the Néel and VBS correlators are
quite small. This is a rather surprising result and does not seem to be due to
any particular small parameter. The magnitude of these anomalous dimensions do
not decrease upon increasing the numerical precision of our integration, so we
believe them to be small but not identically zero. We find that the gauge
fluctuations generally enhance Néel and VBS correlations, whereas Higgs
fluctuations suppress them - the combined result is the stated anomalous
dimensions. As such, we cannot make a strong statement regarding which
ordering the unstable $\operatorname{U}(1)$ phase will prefer, as neither the
Néel nor VBS order parameter show exceptionally enhanced correlations. Higher-
order corrections may show a clearer preference to either Néel or VBS
ordering.

###  IV.4 Tree-level effect of velocity anisotropy on correlation functions

As we have emphasized, one of the key features of this critical theory is that
the emergent Lorentz invariance of the staggered flux phase is broken by the
presence of critical Higgs fields, leading to a non-zero value of the
symmetry-allowed velocity anisotropy term. This anisotropy term also has the
effect of breaking the emergent $\operatorname{SU}(4)$ flavor symmetry. We
refer to Ref. [34] for a more extensive study of the intertwining physical
order parameters of the $\operatorname{SU}(4)$ theory and which relations hold
in the presence of the velocity anisotropy - for our purposes, we note that
the emergent $\operatorname{SO}(5)\subset\operatorname{SU}(4)$ symmetry that
relates the Néel and VBS order parameters is broken down to the microscopic
$\operatorname{SO}(3)\times C_{4}$. At tree-level, the scaling dimensions of
the two order parameters are still the same, but the angular profile of their
correlation functions are modified due to the velocity anisotropy. This lack
of an emergent $\operatorname{SO}(2)$ spatial rotation symmetry in the Néel
and VBS correlation functions may be useful as a numerical probe of the
critical behavior, so we study the angular profile in more detail.

We analytically compute the spatial profile of the Néel order parameter at
tree level. This calculation turns out to be feasible non-perturbatively in
the velocity anisotropy $\Phi$. The VBS correlator is more difficult to study
non-perturbatively in the velocity anisotropy, and we will later compute
corrections to leading order in $\Phi$.

The two-point function in momentum space is given by, with $Q(p)$ the fermion
propagator,

$\displaystyle G_{\text{N\'{e}el}}(k)$
$\displaystyle=-\int\frac{\differential[3]{p}}{(2\pi)^{3}}\Tr\left[Q(p)\mu^{y}Q(p+k)\mu^{y}\right]$
(29)
$\displaystyle=-\frac{2}{1-\Phi^{2}}\sum_{a=\pm}\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{p_{0}(p_{0}+k_{0})+ap_{x}(p_{x}+k_{x,a})+ap_{y}(p_{y}+k_{y,a})}{p^{2}(p+k_{\pm})^{2}}$
$\displaystyle=\frac{1}{8(1-\Phi^{2})}\left(\absolutevalue{k_{+}}+\absolutevalue{k_{-}}\right)\,.$

As before, we define $k_{\pm}=(k_{0},k_{x}\pm\Phi k_{y},k_{y}\pm\Phi k_{x})$.
The Fourier transform

$\displaystyle\int\frac{\differential[3]{k}}{(2\pi)^{3}}e^{ik\cdot
r}\absolutevalue{k_{\pm}}$ (30)

can be performed by a change of variables to give

$\displaystyle G_{\text{N\'{e}el}}(r)$
$\displaystyle\sim\sum_{a=\pm}f_{a}(\theta,\phi)^{-2}\frac{1}{r^{4}}$ (31)

where we change to spherical coordinates, $t=r\cos\theta$,
$x=r\sin\theta\sin\phi$, $y=r\sin\theta\cos\phi$, and

$\displaystyle f_{\pm}(\theta,\phi)$
$\displaystyle=1+\frac{\sin^{2}\theta\left(\pm 2\Phi\sin
2\phi+3\Phi^{2}+\Phi^{4}\right)}{(1-\Phi^{2})^{2}}\,.$ (32)

Therefore,

$\displaystyle G_{\text{N\'{e}el}}(r)\sim g(\theta,\phi)\frac{1}{r^{4}}$ (33)

where
$g(\theta,\phi)={1}/{f_{+}(\theta,\phi)^{2}}+{1}/{f_{-}(\theta,\phi)^{2}}$ is
plotted in Fig. 6. We note the enhanced correlations of the Néel order
parameter along the diagonals, which holds true for generic values of $\Phi$.

Figure 6: Plotted are the angular profiles of the equal-time Néel and VBS
correlation functions in real space, at the fixed point value of velocity
anisotropy $\Phi_{c}$. The Néel order parameter shows enhanced correlations
along the diagonals, whereas the VBS correlations are more enhanced along the
$x$ and $y$ directions. Note that we only plot the absolute value of the
correlation function, and the signs of the perturbative Néel and VBS
correlators flip when moving from the $x$ and $y$ axes to the diagonals. As
this feature is not present in the non-perturbative Néel correlator, it is
possible that this feature similarly vanishes at higher orders for the VBS
correlator.

An analogous calculation of the VBS order parameter is less analytically
tractable, as the one-loop integral cannot be made isotropic by a coordinate
transformation. As such, we resort to a perturbative study of the velocity
anisotropy. This gives

$\displaystyle
G_{\text{VBS}}(k)\sim\absolutevalue{k}-\Phi^{2}\left[2\absolutevalue{k}+\frac{k_{x}^{2}k_{y}^{2}-k_{0}^{2}k^{2}}{\absolutevalue{k}^{3}}\right]+\order{\Phi^{4}}$
(34)

in momentum space, or

$\displaystyle G_{\text{VBS}}(r)\sim\frac{1}{r^{4}}\left[8+\Phi^{2}(\cos
2\theta(40+12\cos 4\phi)+\cos 4\theta(6-3\cos 4\phi)-18\cos 4\phi-14)\right]$
(35)

in real space. The equal-time VBS correlation function is plotted in Fig. 6,
showing enhanced correlations along the cardinal directions. Note that the
correlation function changes sign on the diagonals - this is an unusual
feature, and would seemingly indicate lines in real space where the VBS
correlator vanishes. This feature is also present in the $\order{\Phi^{2}}$
corrections to the Néel correlator but ultimately vanishes in the non-
perturbative result, so this result may only be an artifact of the
perturbative expansion. Further details on this calculation can be found in
Appendix E.

##  V Monopoles

On the square lattice, there exists a monopole operator in the staggered flux
phase - the trivial monopole \- that is invariant under all square lattice
symmetries, and hence is an allowed perturbation. To leading order, the
scaling dimension of the monopole operator scales with the number of fermions
and becomes irrelevant for $N_{f}\geq 3$. Hence, the staggered flux phase by
itself is unstable to monopole proliferation - this is the mechanism which we
claim gives rise to ordered phases in the staggered flux phase, as
condensation of the trivial monopole is conjectured to lead to a fermion
chiral mass generation corresponding to either Néel or VBS order [40, 41]. Our
calculations of the critical theory rely on the assumption that the presence
of massless scalar fields screens monopoles and renders them irrelevant at the
critical point. Here, we draw attention to an additional contribution to the
monopole scaling dimension, which is the non-zero anisotropy in the spinon
dispersion. Prior calculations of monopole scaling dimensions in
$\text{QED}_{3}$ assume a Lorentz-invariant action for the fermions, which is
natural in pure $\text{QED}_{3}$ given that velocity anisotropy terms are
irrelevant in a ${1}/{N_{f}}$ expansion. However, as we have shown, the
presence of critical Higgs fields can give rise to a non-zero velocity
anisotropy at the critical point. An important question is whether this
anisotropy increases or decreases the monopole scaling dimension. In contrast
with the direct modification arising from the critical fields, which is
$\order{1}$, the effect of the anisotropy on the monopole scaling dimension is
$\order{N_{f}}$. Such a modification, if calculated perturbatively in $\Phi$,
arises at $\order{\Phi^{2}}$ \- this is still an appreciable shift given the
relatively large anisotropy $\Phi_{c}\approx 0.46$. Previous works have
studied the effects of a spin Hall mass on the monopole scaling dimension [44,
45], although this perturbation is more tractable as the spin Hall mass is
diagonal in the basis of spinor monopole harmonics. In Appendix F, we outline
the structure of a perturbative calculation for calculating the
$\order{\Phi^{2}}$ corrections to the monopole scaling dimension. An important
observation which makes this calculation tractable is that, while the saddle-
point monopole gauge configuration will not take the form of the rotationally-
invariant Dirac monopole, corrections to the scaling dimension arising from
this difference only arise at higher orders in $\Phi$; hence, to lowest non-
trivial order, one can assume a Dirac monopole background. This calculation
ultimately yields a divergent summation of terms involving Wigner $3$-$j$
symbols; we leave for future work further study of how to properly regularize
this calculation.

Additionally, we briefly comment on the relation between this velocity
anisotropy and the monopole quantum numbers. Prior studies on the effects of a
spin Hall mass [44, 45] have found that the presence of such a term induces a
spin polarization on the monopoles. Each fermion flavor has a zero mode in the
presence of a monopole, and half of these zero modes must be filled in order
to maintain gauge neutrality of the monopole. The presence of a spin Hall mass
polarizes these zero modes, which in turn causes a splitting in the scaling
dimension of the monopoles, with the most-relevant monopole being spin
polarized. One may wonder whether a similar valley polarization can arise due
to our velocity anisotropy term due to the presence of a $\mu^{y}$ in the
anisotropy; however, we check that the first-order energy splitting of the
fermion zero modes due to the velocity anisotropy vanishes. Higher-order
corrections including corrections from Higgs and gauge fluctuations will in
general break the six-fold degeneracy of monopole scaling dimensions; in
particular, the monopoles with Néel and VBS quantum numbers will have
different scaling dimensions, which may cause a preference towards a
particular type of symmetry-breaking in the staggered flux phase. Further
study of the spectrum of monopoles at this critical point may be useful for
determining the IR fate of the proximate staggered flux phase - our
observation is that this behavior is more complicated than a simple valley
polarization of the monopoles.

##  VI Conclusions and future directions

We have presented a large-$N_{f}$ analysis of a deconfined critical theory
separating a gapless $\mathbb{Z}_{2}$ spin liquid Z2Azz13 from the
$\operatorname{U}(1)$ staggered flux phase, the latter of which we assume to
be unstable to either Néel or VBS ordering on the square lattice. This
completes the large-$N$ study of the phase diagram shown in Fig. 2, where the
gapless $\mathbb{Z}_{2}$ spin liquid may emerge as a Higgsed phase from either
a $\operatorname{U}(1)$ or $\operatorname{SU}(2)$ gauge theory. Both of these
parent gauge theories are conjectured to be unstable on the square lattice,
and hence we propose the trajectory through the phase diagram as shown in Fig.
2 as a description of the $J_{1}$-$J_{2}$ square lattice antiferromagnet,
where numerical studies suggest a gapless $\mathbb{Z}_{2}$ spin liquid
emerging between Néel and VBS phases.

Our calculations yield several predictions which may be investigated by future
numerical studies. One of the most striking features of both the
$\operatorname{SU}(2)\rightarrow\mathbb{Z}_{2}$ and
$\operatorname{U}(1)\rightarrow\mathbb{Z}_{2}$ transitions is the lack of
Lorentz invariance, spatial rotation invariance (aside from the discrete
$C_{4}$ rotational symmetry), and in the case of the
$\operatorname{SU}(2)\rightarrow\mathbb{Z}_{2}$ transition, even a lack of
traditional scale invariance. The lack of scale invariance takes the form of
correlation functions decaying as $e^{-\ln^{2}(r)}$ rather than power law, the
difference of which is difficult to detect for small system sizes. As such, it
may be more promising to search for lack of Lorentz invariance ($z\neq 1$), or
lack of a full $\operatorname{SO}(2)$ spatial rotation invariance of
correlation functions. We draw attention to the angular profiles of the Néel
and VBS correlation functions shown in Fig. 6, which come from a mean-field
description of the staggered flux state with the inclusion of a symmetry-
allowed velocity anisotropy and predict enhanced Néel correlations along the
diagonals, and enhanced VBS correlations along the cardinal directions.

##  VII Acknowledgements

We are grateful to Zheng-Cheng Gu, Patrick Ledwith, and Alex Thomson for
insightful discussions. This research was supported by the National Science
Foundation under Grant No. DMR-2002850. This work was also supported by the
Simons Collaboration on Ultra- Quantum Matter, which is a grant from the
Simons Foundation (651440, S.S.).

##  Appendix A Alternative proximate $\mathbb{Z}_{2}$ spin liquids

The phase diagram of the theory we have defined contains three phases on a
mean-field level - the $\operatorname{SU}(2)$ spin liquid of the $\pi$-flux
phase, the $\operatorname{U}(1)$ staggered flux phase, and the gapless
$\mathbb{Z}_{2}$ spin liquid whose projective symmetry group labels it as
Z2A$zz$13 according to Wen’s classification [28]. Although these spin liquids
are the ones we believe to be of relevance to the $J_{1}$-$J_{2}$ model,
additional $\mathbb{Z}_{2}$ spin liquids are accessible in this general
framework by a modification of the Higgs couplings [46]. Although there are
many $\mathbb{Z}_{2}$ spin liquids accessible starting from the $\pi$-flux
phase and we will not attempt to study all of them, we will identify the
continuum theories associated with the eight $\mathbb{Z}_{2}$ spin liquids
proximate to both the $\pi$-flux and staggered flux phase.

###  A.1 Symmetry fractionalization in the continuum staggered flux theory

The $\operatorname{SU}(2)$ gauge symmetry introduced in the fermionic spinon
theory of spin liquids means that the microscopic symmetries of the square
lattice need no longer be realized explicitely, but may be realized
projectively, i.e., up to an overall $\operatorname{SU}(2)$ gauge
transformation. This concept holds true in our continuum theory, and we will
use this to identify possible $\mathbb{Z}_{2}$ spin liquids based off of how
the microscopic symmetries are realized. To see this, consider the field
theory describing the transition from the $\pi$-flux phase to the staggered
flux phase, which consists of four Dirac fermions and a three-component
adjoint Higgs field $\Phi_{3}$, both minimally coupled to an
$\operatorname{SU}(2)$ gauge field, and an additional Yukawa coupling

$\displaystyle\Phi_{3}^{a}\overline{\psi}\sigma^{a}\mu^{y}(\gamma^{x}i\partial_{y}+\gamma^{y}i\partial_{x})\psi\equiv\Phi_{3}^{a}\overline{\psi}\sigma^{a}M\psi$
(36)

When $\Phi_{3}$ is uncondensed, it can be integrated out, and the Yukawa
couplings generate terms irrelevant at long distances. When $\Phi_{3}$ is
condensed - for concreteness, $\langle\Phi_{3}^{a}\rangle=\Phi\delta_{az}$ \-
we replace $\Phi_{3}^{a}$ by its expectation value, and the fermion bilinear
$\overline{\psi}\sigma^{z}\mu^{y}(\gamma^{x}i\partial_{y}+\gamma^{y}i\partial_{x})\psi$
naively breaks translation and rotation symmetry according the symmetry
transformations in Table 1. However, the fermion bilinear is invariant under a
combination of microscopic symmetries and $\operatorname{SU}(2)$ gauge
transformations:

$\displaystyle\mathscr{P}G\,:\,\overline{\psi}\sigma^{a}M\psi\rightarrow\overline{\psi}U_{G}^{\dagger}V_{G}^{\dagger}\sigma^{a}V_{G}MU_{G}\psi\,.$
(37)

$U_{G}$ is the action of the microscopic symmetry and $V_{G}$ is a gauge
transformation, given by

$\displaystyle V_{tx}$ $\displaystyle=g(\phi_{x})i\sigma^{x},$ $\displaystyle
V_{px}$ $\displaystyle=g(\phi_{px}),$ $\displaystyle V_{r}$
$\displaystyle=g(\phi_{r})i\sigma^{x},$ (38) $\displaystyle V_{ty}$
$\displaystyle=g(\phi_{y})i\sigma^{x},$ $\displaystyle V_{py}$
$\displaystyle=g(\phi_{py}),$ $\displaystyle V_{t}$
$\displaystyle=g(\phi_{t}),$

where $g(\phi)\equiv e^{i\phi\sigma_{z}}$ reflects the residual
$\operatorname{U}(1)$ symmetry. The convention that we use for these
subscripts are as follows. Translations along one lattice site in the $x$, $y$
direction are represented by subscripts $tx$, $ty$, respectively. Reflections
along the $x$, $y$ axis are indicated by $px$, $py$. The subscript $r$
corresponds to a $\pi/2$ rotation, and $t$ indicates time-reversal symmetry.
These projective symmetry transformations can be deduced by the way the
microscopic symmetries act on the Dirac fermions, given in Table 1. The
ability to condense additional Higgs fields to yield a $\mathbb{Z}_{2}$ spin
liquid is constrained by the requirement that there exists a choice of phases
$\phi_{i}$ such that all bilinear terms are invariant under the same
projective symmetry transformation. We use this fact to identify the Yukawa
couplings that correspond to the different possible proximate spin liquids -
the phases $\phi_{i}$ can be read off via the symmetry fractionalization of
the spin liquids, as shown in the next section, which uniquely identifies the
Yukawa couplings consistent with these phases.

###  A.2 Translationally-invariant spin liquid ansatzes

There are four proximate $\mathbb{Z}_{2}$ spin liquids with translationally-
invariant mean-field ansatzes - i.e., there exists a gauge in which
translational symmetry is realized explicitly in the mean-field ansatz. These
spin liquids are Z2A$zz$13, Z2A00$1n$, Z2A$zz1n$, and Z2A0013. These spin
liquids are classified based off of their symmetry fractionalization; in other
words, how symmetry operations like $T_{y}^{-1}T_{x}T_{y}T_{x}^{-1}$ are not
directly equivalent to the identity, but equivalent up to a gauge
transformation. These values are shown in Table 2, along with the values
corresponding to the continuum $\operatorname{U}(1)$ staggered flux phase,
which can be read off from the transformations 38.

Table 2: Symmetry fractionalization of the translationally-invariant $\mathbb{Z}_{2}$ spin liquids proximate to the staggered flux phase. | Group Relation | Staggered Flux | Z2A$zz$13 | Z2A0013 | Z2A$zz1n$ | Z2A00$1n$
---|---|---|---|---|---|---
1 | $T_{y}^{-1}T_{x}T_{y}T_{x}^{-1}$ | $-e^{-2i(\phi_{x}-\phi_{y})\sigma^{z}}$ | 1 | 1 | 1 | 1
2 | $P_{y}^{-1}T_{x}P_{y}T_{x}^{-1}$ | $e^{2i\phi_{py}\sigma^{z}}$ | -1 | 1 | -1 | 1
3 | $P_{y}^{-1}T_{y}P_{y}T_{y}$ | $e^{2i\phi_{py}\sigma^{z}}$ | -1 | 1 | -1 | 1
4 | $P_{y}^{2}$ | $e^{2i\phi_{py}\sigma^{z}}$ | -1 | 1 | -1 | 1
5 | $P_{y}^{-1}R_{\pi/2}P_{y}R_{\pi/2}^{-1}$ | $-e^{2i\phi_{py}\sigma^{z}}$ | 1 | -1 | 1 | -1
6 | $R_{\pi/2}^{4}$ | $1$ | 1 | 1 | 1 | 1
7 | $R_{\pi/2}^{-1}T_{x}R_{\pi/2}T_{y}$ | $e^{i(2\phi_{r}-\phi_{x}-\phi_{y})\sigma^{z}}$ | -1 | 1 | -1 | 1
8 | $R_{\pi/2}^{-1}T_{y}R_{\pi/2}T_{x}^{-1}$ | $e^{i(2\phi_{r}-\phi_{x}-\phi_{y})\sigma^{z}}$ | -1 | 1 | -1 | 1
9 | $\mathcal{T}^{-1}R_{\pi/2}^{-1}\mathcal{T}R_{\pi/2}$ | $e^{-2i\phi_{t}\sigma^{z}}$ | -1 | -1 | 1 | 1
10 | $\mathcal{T}^{-1}P_{y}^{-1}\mathcal{T}P_{y}$ | 1 | 1 | 1 | 1 | 1
11 | $\mathcal{T}^{-1}T_{x}^{-1}\mathcal{T}T_{x}$ | $-e^{-2i\phi_{t}\sigma^{z}}$ | 1 | 1 | -1 | -1
12 | $\mathcal{T}^{-1}T_{y}^{-1}\mathcal{T}T_{y}$ | $-e^{-2i\phi_{t}\sigma^{z}}$ | 1 | 1 | -1 | -1
13 | $\mathcal{T}^{2}$ | $e^{-2i\phi_{t}\sigma^{z}}$ | -1 | -1 | 1 | 1

Using this table to identify the phases $\phi_{i}$ we find the continuum
projective symmetry group for Z2A$zz$13,

$\displaystyle V_{tx}$ $\displaystyle=-i\sigma^{y}\,,\quad V_{px}=\pm
i\sigma^{z}$ (39) $\displaystyle V_{r}$
$\displaystyle=-\frac{i}{\sqrt{2}}(\sigma^{x}-\sigma^{y})$ $\displaystyle
V_{ty}$ $\displaystyle=-i\sigma^{x}\,,\quad V_{py}=-i\sigma^{z}\,,\quad
V_{t}=i\sigma^{z}$

Z2Azz1n,

$\displaystyle V_{tx}$ $\displaystyle=-i\sigma^{y}\,,\quad V_{px}=\pm
i\sigma^{z}$ (40) $\displaystyle V_{r}$
$\displaystyle=-\frac{i}{\sqrt{2}}(\sigma^{x}-\sigma^{y})$ $\displaystyle
V_{ty}$ $\displaystyle=-i\sigma^{x}\,,\quad V_{py}=-i\sigma^{z}\,,\quad
V_{t}=1$

Z2A0013,

$\displaystyle V_{tx}$ $\displaystyle=-i\sigma^{y}\,,\quad V_{px}=1$ (41)
$\displaystyle V_{r}$
$\displaystyle=-\frac{i}{\sqrt{2}}(\sigma^{x}+\sigma^{y})$ $\displaystyle
V_{ty}$ $\displaystyle=-i\sigma^{x}\,,\quad V_{py}=1\,,\quad
V_{t}=i\sigma^{z}$

and Z2A001n,

$\displaystyle V_{tx}$ $\displaystyle=-i\sigma^{y}\,,\quad V_{px}=1$ (42)
$\displaystyle V_{r}$
$\displaystyle=-\frac{i}{\sqrt{2}}(\sigma^{x}+\sigma^{y})$ $\displaystyle
V_{ty}$ $\displaystyle=-i\sigma^{x}\,,\quad V_{py}=1\,,\quad V_{t}=1\,.$

From this, we can identify the Yukawa couplings consistent with these
transformations, given in Table 3. In order to realize these transitions, we
require two Higgs fields, $\Phi_{1}$ and $\Phi_{2}$ which transform into each
other under square lattice rotations and condense to acquire a non-zero
$\langle\Phi_{1}^{x}\rangle$ and $\langle\Phi_{2}^{y}\rangle$. The fact that
these fields are related under the microscopic rotational symmetry requires
both Higgs fields to have the same mass, and hence they can both be condensed
simultaneously by tuning a single parameter. These theories may be studied in
a manner analogous to our study of the Z2A$zz$13 transition.

Table 3: Yukawa couplings of the two adjoint Higgs fields which realize the symmetry fractionalization of the proximate translationally-invariant spin liquids. | Z2A$zz$13 | Z2A0013 | Z2A$zz1n$ | Z2A00$1n$
---|---|---|---|---
$\Phi_{1}$ coupling | $\mu^{x}\gamma^{y}$ | $\mu^{x}\gamma^{x}$ | $\mu^{x}i\partial_{x}$ | $\mu^{z}i\partial_{x}$
$\Phi_{2}$ coupling | $\mu^{z}\gamma^{x}$ | $\mu^{z}\gamma^{y}$ | $\mu^{z}i\partial_{y}$ | $\mu^{x}i\partial_{y}$

###  A.3 Non-translationally-invariant spin liquid ansatzes

The four additional spin liquid phases proximate to the staggered flux phase
are Z2Bzz13, Z2B0013, Z2B001n, and Z2Bzz1n. These spin are distinct from the
first four as one cannot write down a translationally-invariant mean-field
ansatz for them - note that this does not correspond to a physical breaking of
the translational symmetry, as the symmetry is still realized projectively;
rather, the statement is that the symmetry must be realized projectively. The
symmetry fractionalization of these spin liquids is identical to their Z2A
counterparts in Table 2 aside from a change in sign in rows 1, 7, and 8. This
leads to the continuum projective symmetry group for Z2Bzz13,

$\displaystyle V_{tx}$ $\displaystyle=i\sigma^{x}\,,\quad V_{px}=\pm
i\sigma^{z}$ (43) $\displaystyle V_{r}$ $\displaystyle=i\sigma^{x}$
$\displaystyle V_{ty}$ $\displaystyle=i\sigma^{x}\,,\quad
V_{py}=-i\sigma^{z}\,,\quad V_{t}=i\sigma^{z}$

Z2Bzz1n,

$\displaystyle V_{tx}$ $\displaystyle=i\sigma^{x}\,,\quad V_{px}=\pm
i\sigma^{z}$ (44) $\displaystyle V_{r}$ $\displaystyle=i\sigma^{x}$
$\displaystyle V_{ty}$ $\displaystyle=i\sigma^{x}\,,\quad
V_{py}=-i\sigma^{z}\,,\quad V_{t}=1$

Z2B0013,

$\displaystyle V_{tx}$ $\displaystyle=i\sigma^{x}\,,\quad V_{px}=1$ (45)
$\displaystyle V_{r}$ $\displaystyle=i\sigma^{y}$ $\displaystyle V_{ty}$
$\displaystyle=i\sigma^{x}\,,\quad V_{py}=1\,,\quad V_{t}=i\sigma^{z}$

and Z2B001n,

$\displaystyle V_{tx}$ $\displaystyle=i\sigma^{x}\,,\quad V_{px}=1$ (46)
$\displaystyle V_{r}$ $\displaystyle=i\sigma^{y}$ $\displaystyle V_{ty}$
$\displaystyle=i\sigma^{x}\,,\quad V_{py}=1\,,\quad V_{t}=1\,.$

Describing these phases as a Higgsed phase of the $\pi$-flux spin liquid turns
out to differ significantly from their Z2A counterparts. The transition from
the staggered flux phase to these spin liquids is driven by the condensation
of a single Higgs field $\Phi_{1}$, which acquires a non-zero expectation
value $\langle\Phi_{1}^{x}\rangle$ in the $\mathbb{Z}_{2}$ phase. This
prevents a direct transition from the $\operatorname{SU}(2)$ $\pi$-flux state
to the $\mathbb{Z}_{2}$ spin liquid, as the masses of two Higgs fields
$\Phi_{1}$ and $\Phi_{3}$ aren’t constrained to be equal by the microscopic
symmetries.

Table 4: Yukawa couplings of the single adjoint Higgs fields which realize the symmetry fractionalization of the proximate Z2B spin liquids. | Z2Bzz13 | Z2B0013 | Z2Bzz1n | Z2B001n
---|---|---|---|---
$\Phi_{1}$ coupling | $\mu^{y}\partial_{x}\partial_{y}(\partial_{x}^{2}+\partial_{y}^{2})$ | $\mu^{y}$ | $i\partial_{0}$ | $\mu^{y}(\gamma^{y}i\partial_{x}-\gamma^{x}i\partial_{y})$

##  Appendix B Emergent subsystem symmetries in the
$\operatorname{SU}(2)\rightarrow\mathbb{Z}_{2}$ transition

In this appendix, we provide a more detailed discussion of the emergent
subsystem symmetry in the $\operatorname{SU}(2)\rightarrow\mathbb{Z}_{2}$
transition present in our theory. This transition is driven by the
simultaneous condensation of the two Higgs fields $\Phi_{1,2}$, whose masses
are fixed to be equal by the microscopic $C_{4}$ rotation symmetry of the
square lattice. The Lagrangian that describes this transition is

$\displaystyle\mathcal{L}$
$\displaystyle=\mathcal{L}_{\psi}+\mathcal{L}_{\Phi}+\mathcal{L}_{\Phi\psi}$
(47) $\displaystyle\mathcal{L}_{\psi}$
$\displaystyle=i\bar{\psi}\gamma^{\mu}\left(\partial_{\mu}-iA^{a}_{\mu}\sigma^{a}\right)\psi\,.$
$\displaystyle\mathcal{L}_{\Phi}$
$\displaystyle=\sum_{i=1}^{3}D_{\mu}\Phi_{i}^{a}D^{\mu}\Phi_{i}^{a}+V(\Phi)$
$\displaystyle\mathcal{L}_{\Phi\psi}$
$\displaystyle=\Phi_{1}^{a}\,\bar{\psi}\mu^{z}\gamma^{x}\sigma^{a}\psi+\Phi_{2}^{a}\,\bar{\psi}\mu^{x}\gamma^{y}\sigma^{a}\psi$

$V(\Phi)$ contains various symmetry-allowed potential terms for $\Phi_{1,2}$,
all of which are irrelevant other than the mass term which we tune to zero at
criticality. Examining this action, we see that the quadratic fermion action
along with the $\Phi_{1}$ Yukawa coupling in invariant under the
transformation

$\displaystyle\psi$ $\displaystyle\rightarrow
e^{if_{a}(x)\mu^{z}\sigma^{a}}\psi$ (48) $\displaystyle\Phi_{1}^{a}$
$\displaystyle\rightarrow U^{-1}_{ab}(x)\Phi_{1}^{b}+\partial_{x}f_{a}(x)$

and likewise for $\Phi_{2}$,

$\displaystyle\psi$ $\displaystyle\rightarrow
e^{ig_{a}(y)\mu^{z}\sigma^{a}}\psi$ (49) $\displaystyle\Phi_{1}^{a}$
$\displaystyle\rightarrow\tilde{U}^{-1}_{ab}(y)\Phi_{1}^{b}+\partial_{y}g_{a}(y)$

with $U(x)$ the $\operatorname{SO}(3)$ rotation corresponding to the adjoint
$\operatorname{SU}(2)$ action of $e^{if_{a}(x)\sigma^{a}}$, i.e.

$U(x)_{ab}\equiv\frac{1}{2}\Tr\left[e^{-if_{c}(x)\sigma^{c}}\sigma^{a}e^{if_{d}(x)\sigma^{d}}\sigma^{b}\right]$
(50)

and likewise for $\tilde{U}(y)$. The simplest way to see the existence of
these symmetries is to consider $\Phi_{1}$ and $\Phi_{2}$ as the $x$ and $y$
components, respectively, of fictitious gauge fields corresponding to the
$\operatorname{SU}(2)$ symmetries $\psi\rightarrow
e^{i\sigma^{a}\mu^{z,x}}\psi$. The Yukawa couplings present in our theory
couple these Higgs fields to the $x$ and $y$ components of the respective
conserved currents. As such, the $\Phi_{1}$ ($\Phi_{2}$) Yukawa coupling and
the fermion quadratic term are invariant under $x$-dependent ($y$-dependent)
$\operatorname{SU}(2)$ transformations in a manner analogous to gauge
invariance. Of course, this symmetry is ultimately broken in the full
Lagrangian, both by the actual gauge field $A_{\mu}$ as well as the Higgs
potential $V(\Phi)$; however, in a $1/N_{f}$ expansion, the leading-order
effective action for the Higgs fields does possess this symmetry. This
effective action is more relevant at long distances than subleading
corrections, such as the bare Higgs action
$\Phi_{1,2}^{a}\partial^{2}\Phi_{1,2}^{a}$, and hence our analysis suggests
that these emergent subsystem symmetries control the behavior of physical
observables at criticality. The $\order{N_{f}^{-1}}$ contributions to the Néel
and VBS order parameters have logarithm squared divergences, which we
conjecture should lead to correlations decaying as
$\frac{1}{r^{\alpha}}e^{-\beta\ln^{2}(r)}$, with
$\beta=-\frac{6}{\pi^{2}N_{f}}$ for the VBS correlator and
$-\frac{12}{\pi^{2}N_{f}}$ for the Néel correlator, and $\alpha$ being some
non-universal coefficient. Interestingly, we note that the logarithm squared
divergences for correlations of the scalar spin chirality exactly cancel at
$\order{N_{f}^{-1}}$. We expect that this is related to the invariance of the
scalar spin chirality order parameter under the subsystem symmetries given by
Eqs. (48) and (49). Many open questions remain in regards to these emergent
symmetries; in particular, the physical interpretation of them is not clear -
while these symmetry transformations are not purely gauge transformations,
they do contain an action on the gauge $\operatorname{SU}(2)$ space given by
the $\sigma^{a}$ operators, and as such, can rotate gauge-invariant operators
such as $\bar{\psi}\mu^{i}\psi$ into non-gauge-invariant ones.

##  Appendix C Calculation of effective bosonic propagators

In order to study our critical theory in a ${1}/{N_{f}}$ expansion, we must
calculate the effective propagators for both the Higgs and gauge bosons, which
are generated by one-loop fermion diagrams. These calculations are complicated
substantially by the non-Lorentz-invariant nature of the fermion propagator.
As given in Eq. (9), the bare fermion action is given by

$\bar{\psi}\left[\not{k}+\Phi\mu^{y}\sigma^{z}(\gamma^{y}k_{x}+\gamma^{x}k_{y})\right]\psi\,.$
(51)

It is useful to work in an eigenbasis of $\sigma^{z}\mu^{y}$. Inverting this,
we get a diagonal $4\times 4$ matrix in gauge/valley space, with elements

$\displaystyle\frac{\gamma^{0}k_{0}+\gamma^{x}(k_{x}\pm\Phi
k_{y})+\gamma^{y}(k_{y}\pm\Phi k_{x})}{k_{0}^{2}+(k_{x}\pm\Phi
k_{y})^{2}+(k_{y}\pm\Phi k_{x})^{2}}\,.$ (52)

As in the main text, we define the variables

$\displaystyle k_{x,\pm}\equiv k_{x}\pm\Phi k_{y}\,,\quad k_{y,\pm}\equiv
k_{y}\pm\Phi
k_{x}\,,\absolutevalue{k_{\pm}}\equiv\sqrt{k_{0}^{2}+k_{x,\pm}^{2}+k_{y,\pm}^{2}}\,.$
(53)

The inverse Higgs propagator is generated by the one-loop diagram given by
Fig. 3 of the main text, which translates to the integral

$\displaystyle-\frac{4}{N_{f}}\sum_{a=\pm}\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{p_{0}(p_{0}+k_{0})+ap_{x,a}(p_{y,a}+k_{y,a})+ap_{y,a}(p_{x,a}+k_{x,a}))}{k_{a}^{2}(k+p)_{a}^{2}}$
(54)

To simplify this integral, we perform a change of integration variables to
$(p_{0},p_{x,\pm},p_{y,\pm})$. Performing this change of variables gives a
factor of $(1-\Phi^{2})^{-1}$ in the integral.

To evaluate this integral, we calculate the two general forms of integrals
relevant to Eq. (54).

$\displaystyle\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{p_{0}(p_{0}+k_{0})}{p^{2}(p+k)^{2}}$
(55)
$\displaystyle=\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{p_{0}(p_{0}+k_{0})}{\left[(p+xk)^{2}+\Delta\right]^{2}}\quad\Delta\equiv
k^{2}x(1-x)$
$\displaystyle=\int_{0}^{1}\differential{x}\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{(p_{0}-xk_{0})(p_{0}+(1-x)k_{0})}{(p^{2}+\Delta)^{2}}$
$\displaystyle=\frac{1}{2\pi^{2}}\int_{0}^{1}\differential{x}\int_{0}^{\infty}\differential{p}\frac{p^{2}[p^{2}/3-x(1-x)k_{0}^{2}]}{(p^{2}+\Delta)^{2}}$
$\displaystyle=\frac{1}{2\pi^{2}}\int_{0}^{1}\differential{x}\left[-\frac{\pi\sqrt{\Delta}}{4}-x(1-x)\frac{k_{0}^{2}\pi}{4\sqrt{\Delta}}\right]$
$\displaystyle=-\frac{k^{2}+k_{0}^{2}}{64\absolutevalue{k}}$
$\displaystyle\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{p_{x}(p_{y}+k_{y})}{p^{2}(p+k)^{2}}$
$\displaystyle=\int_{0}^{1}\differential{x}\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{(p_{x}-xk_{x})(p_{y}+(1-x)k_{y})}{(p^{2}+\Delta)^{2}}$
$\displaystyle=-\int_{0}^{1}\differential{x}\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{x(1-x)k_{x}k_{y}}{(p^{2}+\Delta)^{2}}$
$\displaystyle=-\frac{k_{x}k_{y}}{2\pi^{2}}\int_{0}^{1}\differential{x}\int_{0}^{\infty}\differential{p}\frac{x(1-x)}{(p^{2}+\Delta)^{2}}$
$\displaystyle=-\frac{k_{x}k_{y}}{64\absolutevalue{k}}$

In these calculations, we omit divergent terms which renormalize the boson
mass, as this term is set to zero at criticality. Substituting this expression
back into Eq. (12) for general $p_{i},k_{i}$, we get the Higgs propagator

$\displaystyle\Gamma$
$\displaystyle=\frac{\lambda^{2}}{16N_{f}(1-\Phi^{2})}\left[\frac{k_{+}^{2}+k_{0}^{2}+2k_{x,+}k_{y,+}}{\absolutevalue{k_{+}}}+\frac{k_{-}^{2}+k_{0}^{2}-2k_{x,-}k_{y,-}}{\absolutevalue{k_{-}}}\right]$
(56)

##  Appendix D Derivation of one-loop renormalization group equations

In this appendix, we give a derivation of the renormalization group equations
used in the main text. The one-loop contributions to the fermion self-energy
$\Sigma(k)$ are UV divergent, and hence require a UV cutoff $\Lambda$. The
behavior of the self-energy upon integrating out high-energy modes is dictated
by the logarithmic derivative with respect to the cutoff,
$\Lambda\derivative{\Lambda}\Sigma(k)$. The fact that our propagators are
homogeneous functions of the three-momenta allow us to calculate this
logarithmic derivative explicitly without reference to a specific cutoff. We
assume that our regularized one-loop expression for the self-energy takes the
form

$\displaystyle\Sigma(k)=\int\frac{\differential[3]{p}}{(2\pi)^{3}}F(p+k)G(p)\mathcal{K}\left(\frac{p^{2}}{\Lambda^{2}}\right)\mathcal{K}\left(\frac{(k+p)^{2}}{\Lambda^{2}}\right)$
(57)

where $F$ and $G$ are homogeneous functions of the three-momenta with degree
$-1$ \- we take $F$ to be the fermion propagator, and $G$ to be the boson
propagator (either Higgs or gauge) along with the various vertex coefficients.
The function $\mathcal{K}(y)$ serves as a UV cutoff with the property that
$\mathcal{K}(0)=1$ and $\mathcal{K}(y)$ falls off rapidly for large $y$, i.e.,
$\mathcal{K}(y)=e^{-y}$. Since we are interested in the behavior at small
momenta, we expand around $k=0$,

$\displaystyle\Sigma(k)\approx
k_{\mu}\int\frac{\differential[3]{p}}{(2\pi)^{3}}\left[\partialderivative{F(p)}{p_{\mu}}G(p)\mathcal{K}^{2}\left(\frac{p^{2}}{\Lambda^{2}}\right)+F(p)G(p)\mathcal{K}\left(\frac{p^{2}}{\Lambda^{2}}\right)\frac{2p_{\mu}}{\Lambda^{2}}\mathcal{K}^{\prime}\left(\frac{p^{2}}{\Lambda^{2}}\right)\right]\,.$
(58)

We then take the logarithmic derivative,

$\displaystyle\Lambda\derivative{\Lambda}\Sigma(k)$ $\displaystyle\approx
k_{\mu}\int\frac{\differential[3]{p}}{(2\pi)^{3}}\Bigg{[}\left\\{-\frac{4p^{2}}{\Lambda^{2}}\partialderivative{F(p)}{p_{\mu}}-4F(p)\frac{p_{\mu}}{\Lambda^{2}}\right\\}G(p)\mathcal{K}\left(\frac{p^{2}}{\Lambda^{2}}\right)\mathcal{K}^{\prime}\left(\frac{p^{2}}{\Lambda^{2}}\right)$
(59)
$\displaystyle-4\frac{p^{2}p_{\mu}}{\Lambda^{4}}F(p)G(p)\left\\{\mathcal{K}\left(\frac{p^{2}}{\Lambda^{2}}\right)\mathcal{K}^{\prime\prime}\left(\frac{p^{2}}{\Lambda^{2}}\right)+\mathcal{K}^{\prime
2}\left(\frac{p^{2}}{\Lambda^{2}}\right)\right\\}\Bigg{]}\,.$

We now convert to spherical coordinates,
$\vec{p}=y\Lambda(\cos\theta,\sin\theta\sin\phi,\sin\theta\cos\phi)$, and use
the homogeneity property of $F$ and $G$ to pull out factors of
$(y\Lambda)^{-1}$.

$\displaystyle\Lambda\frac{d}{d\Lambda}\Sigma(k)$ $\displaystyle\approx$
$\displaystyle\frac{k_{\mu}}{8\pi^{3}}\int_{0}^{\pi}\sin\theta\differential{\theta}\int_{0}^{2\pi}\differential{\phi}\left[\left\\{-4\frac{\partial
F(\hat{p})}{\partial
p_{\mu}}-4\hat{p}_{\mu}F(\hat{p})\right\\}G(\hat{p})\int_{0}^{\infty}y\differential{y}\mathcal{K}(y^{2})\mathcal{K}^{\prime}(y^{2})\right.$
(60)
$\displaystyle\left.~{}~{}-4\hat{p}_{\mu}F(\hat{p})G(\hat{p})\int_{0}^{\infty}y^{3}dy\left\\{\mathcal{K}\left(y^{2}\right)\mathcal{K}^{\prime\prime}\left(y^{2}\right)+\mathcal{K}^{\prime
2}\left(y^{2}\right)\right\\}\right]$

The integral over $y$ can be done explicitly via integration by parts, which
causes the dependence on the cutoff function $\mathcal{K}$ to drop out. This
leads to the expression cited in the main text

$\displaystyle\Lambda\derivative{\Lambda}\Sigma(k)$
$\displaystyle=\frac{k_{\lambda}}{8\pi^{3}}\int_{0}^{2\pi}\differential{\phi}\int_{0}^{\pi}\sin\theta\differential{\theta}\partialderivative{F(\hat{p})}{p_{\lambda}}G(\hat{p})\,.$
(61)

where $\hat{p}\equiv(\cos\theta,\sin\theta\sin\phi,\sin\theta\cos\phi)$.

Explicitly, we take, for the Higgs contribution to the self-energy, defining
$Q(p)$ as the fermion propagator,

$\displaystyle F(p)$
$\displaystyle=\left(\mu^{z}\gamma^{x}+i\mu^{x}\gamma^{y}\right)\sigma^{-}Q(p)\left(\mu^{z}\gamma^{x}-i\mu^{x}\gamma^{y}\right)\sigma^{+}+\left(\mu^{z}\gamma^{x}-i\mu^{x}\gamma^{y}\right)\sigma^{+}Q(p)\left(\mu^{z}\gamma^{x}+i\mu^{x}\gamma^{y}\right)\sigma^{-}\,,$
(62) $\displaystyle G(p)$ $\displaystyle=\frac{1}{\Gamma(p)}\,.$

For the gauge contribution,

$\displaystyle F_{\mu\nu}(p)$
$\displaystyle=\left[\gamma^{\mu}+\Phi\mu^{y}\sigma^{z}\left(\delta_{\mu,x}\gamma^{y}+\delta_{\mu,y}\gamma^{x}\right)\right]Q(p)\left[\gamma^{\nu}+\Phi\mu^{y}\sigma^{z}\left(\delta_{\nu,x}\gamma^{y}+\delta_{\nu,y}\gamma^{x}\right)\right]\,,$
(63) $\displaystyle G^{\mu\nu}(p)$
$\displaystyle=\left(\Pi^{-1}(p)\right)^{\mu\nu}\,,$

where the fermion self-energy now contains a summation over $\mu\,,\nu$. A
similar approach can be used to regulate the one-loop vertex corrections,
which take the form

$\displaystyle\Xi_{i}=\int\frac{\differential[3]{p}}{(2\pi)^{3}}H_{i}(p)\mathcal{K}^{3}\left(\frac{p^{2}}{\Lambda^{2}}\right)$
(64)

where $H_{i}(p)$ is a homogeneous function of degree $-3$. Upon taking the
logarithmic derivative,

$\displaystyle\Lambda\derivative{\Lambda}\Xi_{i}=-3\int\frac{\differential[3]{p}}{(2\pi)^{3}}H_{i}(p)p^{2}\mathcal{K}^{2}\left(\frac{p^{2}}{\Lambda^{2}}\right)\mathcal{K}^{\prime}\left(\frac{p^{2}}{\Lambda^{2}}\right)$
(65)

Converting to spherical coordinates and integrating by parts, we get

$\Lambda\derivative{\Lambda}\Xi_{i}=\frac{1}{8\pi^{3}}\int_{0}^{2\pi}\differential{\phi}\int_{0}^{\pi}\sin\theta\differential{\theta}H_{i}(\hat{p})\,.$
(66)

Explicitly, the Higgs correction is

$\displaystyle H_{i}(p)$
$\displaystyle=Q(p)\mu^{i}Q(p)\left(\mu^{z}\gamma^{x}+i\mu^{x}\gamma^{y}\right)\sigma^{-}\frac{1}{\Gamma(p)}\left(\mu^{z}\gamma^{x}-i\mu^{x}\gamma^{y}\right)\sigma^{+}$
(67)
$\displaystyle+Q(p)\mu^{i}Q(p)\left(\mu^{z}\gamma^{x}-i\mu^{x}\gamma^{y}\right)\sigma^{+}\frac{1}{\Gamma(p)}\left(\mu^{z}\gamma^{x}+i\mu^{x}\gamma^{y}\right)\sigma^{-}$

and the gauge correction,

$\displaystyle H_{i}(p)$
$\displaystyle=Q(p)\mu^{i}Q(p)\left[\gamma^{\mu}+\Phi\mu^{y}\sigma^{z}\left(\delta_{\mu,x}\gamma^{y}+\delta_{\mu,y}\gamma^{x}\right)\right]\left(\Pi(p)^{-1}\right)^{\mu\nu}\left[\gamma^{\nu}+\Phi\mu^{y}\sigma^{z}\left(\delta_{\nu,x}\gamma^{y}+\delta_{\nu,y}\gamma^{x}\right)\right]$
(68)

##  Appendix E Anisotropic correlation functions in real space

As shown in the main text, the momentum-space Néel correlator is given by

$\displaystyle
G_{\text{N\'{e}el}}(k)=\frac{1}{8(1-\Phi^{2})}\left(\absolutevalue{k_{+}}+\absolutevalue{k_{-}}\right)$
(69)

where we define $\absolutevalue{k_{\pm}}=\sqrt{k_{0}^{2}+(k_{x}\pm\Phi
k_{y})^{2}+(k_{y}\pm\Phi k_{x})^{2}}$. The Fourier transform of this function
can be computed with the knowledge of the (suitably regularized) Fourier
transform in three dimensions, $\absolutevalue{k}\rightarrow\frac{1}{r^{4}}$.
We take the Fourier transform

$\displaystyle\int\frac{\differential[3]{k}}{(2\pi)^{3}}e^{ik\cdot
r}\absolutevalue{k_{\pm}}$ (70)

and perform a change of variables to shift the anisotropy to the spatial
coordinates

$\displaystyle\frac{1}{1-\Phi^{2}}\int\frac{\differential[3]{k}}{(2\pi)^{3}}e^{ik\cdot
r_{\pm}}\absolutevalue{k}$ (71) $\displaystyle t_{\pm}$ $\displaystyle=t$
$\displaystyle x_{\pm}$ $\displaystyle=\frac{1}{1-\Phi^{2}}(x\pm\Phi y)$
$\displaystyle y_{\pm}$ $\displaystyle=\frac{1}{1-\Phi^{2}}(y\pm\Phi x)$

which yields the real space correlator given in the main text.

To compute the Fourier transform of the VBS correlator, given perturbatively
by

$\displaystyle
G_{\text{VBS}}(k)\sim\absolutevalue{k}-\Phi^{2}\left[2\absolutevalue{k}+\frac{k_{x}^{2}k_{y}^{2}-k_{0}^{2}k^{2}}{\absolutevalue{k}^{3}}\right]+\order{\Phi^{4}}\,,$
(72)

we define the function

$\displaystyle
f(a_{i},k_{i})=\sqrt{a_{0}k_{0}^{2}+a_{x}k_{x}^{2}+a_{y}k_{y}^{2}}\,.$ (73)

The Fourier transform of this function can be calculated by a similar change
of variables,

$\displaystyle\int\frac{\differential[3]{k}}{(2\pi)^{3}}f(a_{i},k_{i})\sim\frac{1}{\sqrt{a_{0}a_{x}a_{y}}}\left(\frac{t^{2}}{a_{0}}+\frac{x^{2}}{a_{x}}+\frac{y^{2}}{a_{y}}\right)^{-2}\,.$
(74)

The various terms in the $\order{\Phi^{2}}$ corrections to the VBS correlator
can be obtained by taking derivatives of $f(a_{i},k_{i})$ with respect to
$a_{i}$ and setting $a_{i}=1$. This allows us to calculate the real space VBS
correlator and gives the result in the main text.

##  Appendix F Perturbative corrections to monopole scaling dimension

In this appendix, we present a partial calculation of the $\order{N_{f}}$
corrections to the scaling dimensions of a monopole at our deconfined critical
point. As previously established [37], the $\order{N_{f}}$ scaling dimension
for isotropic $\text{QED}_{3}$ is $\Delta=1.06N_{f}$. We present the scaling
dimension using our convention, where $\text{QED}_{3}$ with $N=4$ Dirac
fermions corresponds to $N_{f}=1$. Although gauge and Higgs fluctuations give
corrections to this value, these corrections are subleading in $N_{f}$, and
the only $\order{N_{f}}$ correction comes from taking the saddle-point
solutions of the bosonic fields and calculating the shift in free energy
arising from the anisotropic Dirac dispersion relation. We proceed
perturbatively in the Dirac anisotropy parameter $\Phi$ \- this is necessary
as the anisotropy will in principle modify the saddle-point monopole
configuration of the gauge field. As we will see, to leading order in $\Phi$,
the gauge field configuration corresponding to the isotropic Dirac monopole
will be sufficient.

We start with the action for $\text{QED}_{3}$ with the allowed velocity
anisotropy term, omitting the Higgs fields as they will not play any role in
the calculation

$S=i\int\differential[3]{r}\left[\bar{\psi}\not{D}\psi+\Phi\bar{\psi}\mu^{y}\left(\gamma^{x}D_{x}-\gamma^{y}D_{y}\right)\right]\psi\,.$
(75)

We leave implicit the summation over the $4N_{f}$ fermions. Note that this
action is different than in the main text. This is because we follow the
convention used in [34], where the gauge field is coupled in the usual way,
$D_{\mu}\equiv\partial_{\mu}-iA_{\mu}$, and the microscopic
$\operatorname{SU}(2)$ spin rotation symmetry is implemented explicitly by the
$\sigma^{i}$ matrices. We refrain from using this convention in the main
calculation, as the coupling to the Higgs field is not easily expressible in
this form and overall makes the calculation more complicated.

In the absence of a velocity anisotropy, the saddle-point configurations for
the gauge field corresponding to $n$ units of magnetic flux at the origin are
given by

$\bar{A}_{n}(r)=\frac{n}{2}(1-\cos\theta)\differential{\phi}$ (76)

The non-zero anisotropy will affect these saddle-point solutions. The leading
order corrections to these solutions are $\order{\Phi^{2}}$, as the
$\order{N_{f}}$ effective action for the gauge field upon integrating out the
fermions has only corrections of $\order{\Phi^{2}}$ and higher. Hence, we
write the saddle-point gauge field solution in the presence of a velocity
anisotropy as $A_{n}(r)\equiv\bar{A}_{n}(r)+\delta A_{n}$, with $\delta
A_{n}\sim\order{\Phi^{2}}$.

In order to calculate the scaling dimension of this monopole, we set
$r=e^{\tau}$ and perform a Weyl rescaling

$\displaystyle g_{\mu\nu}$ $\displaystyle\rightarrow e^{-2\tau}g_{\mu\nu}$
(77) $\displaystyle\bar{\psi}\,,\psi$ $\displaystyle\rightarrow
e^{-\tau}\bar{\psi}\,,e^{-\tau}\psi$

This rescaling maps the scaling dimension of the monopole operator to the free
energy $F=-\log Z$ of the system [47].

To leading order in $N_{f}$, we ignore gauge and Higgs fluctuations, and the
action reduces down to one of free fermions with a background gauge field.
This action can be put in a nearly-diagonal form with the aid of monopole
harmonics [48] and their spinor generalization [37]. By expanding $\psi$ in
terms of these harmonics,

$\psi(r)=\int\frac{\differential{\omega}}{2\pi}\left[\sum_{\ell=n/2}^{\infty}\sum_{m=-\ell-1}^{\ell}\Psi_{T}^{\ell
m}(\omega)T_{n,\ell
m}(\theta,\phi)+\sum_{\ell=n/2}^{\infty}\sum_{m=-\ell}^{\ell}\Psi_{S}^{\ell
m}(\omega)S_{n,\ell m}(\theta,\phi)\right]e^{-i\omega\tau}$ (78)

where $T_{n,\ell m}\,,S_{n,\ell m}$ are eigenvalues of the orbital angular
momentum operator $\vec{L}^{2}$ in the presence of a strength $n$ monopole,
with orbital angular momentum $\ell$ and total angular momentum $\ell+1/2$ for
$T_{n,\ell m}$ and $\ell-1/2$ for $S_{n,\ell m}$. Explicit expressions for the
spinor harmonics $T_{n,\ell m}$ and $S_{n,\ell m}$ may be found in [37]. The
variables $\Psi_{T}^{\ell m}\,,\Psi_{S}^{\ell m}$ are anti-commuting
coefficients.

Expanded in this form, the isotropic action with $\Phi=0$ can be written as

$S_{0}=\int\frac{\differential{\omega}}{2\pi}\sum_{\ell=n/2}^{\infty}\sum_{m=-\ell}^{\ell-1}\begin{pmatrix}\Psi_{T}^{(\ell-1)m}(\omega)^{*}&\Psi_{S}^{\ell
m}(\omega)^{*}\end{pmatrix}\mathbf{N}_{n,\ell}(\omega+i\mathbf{M}_{n,\ell})\begin{pmatrix}\Psi_{T}^{(\ell-1)m}(\omega)\\\
\Psi_{S}^{\ell m}(\omega)\end{pmatrix}$ (79)

with

$\displaystyle\mathbf{M}_{n,\ell}$
$\displaystyle=\begin{pmatrix}\ell\left(1-\frac{n^{2}}{4\ell^{2}}\right)&-\frac{n}{2}\sqrt{1-\frac{n^{2}}{4\ell^{2}}}\\\
-\frac{n}{2}\sqrt{1-\frac{n^{2}}{4\ell^{2}}}&\ell\left(1-\frac{n^{2}}{4\ell^{2}}\right)\end{pmatrix}\,,$
(80) $\displaystyle\mathbf{N}_{n,\ell}$
$\displaystyle=\begin{pmatrix}-\frac{n}{2\ell}&-\sqrt{1-\frac{n^{2}}{4\ell^{2}}}\\\
-\sqrt{1-\frac{n^{2}}{4\ell^{2}}}&\frac{n}{2\ell}\end{pmatrix}\,.$

The free energy is then given by

$\log
Z_{0}=4N_{f}\int\frac{\differential{\omega}}{2\pi}\sum_{\ell=n/2}^{\infty}\log\det(\mathbf{N}_{n,\ell}(\omega+i\mathbf{M}_{n,\ell}))\,.$
(81)

This expression can be evaluated via zeta function regularization, yielding
the aforementioned scaling dimension of $\Delta=1.06N_{f}$.

Corrections to the free energy as a consequence of a non-zero $\Phi$ can be
calculated perturbatively. Writing the action as $S=S_{0}+\delta S$, with
$\delta S\sim\order{\Phi}$, we have

$\displaystyle\log Z$
$\displaystyle=\log\left[\int\mathcal{D}\psi\mathcal{D}\bar{\psi}e^{-S}\right]=\log\left[Z_{0}+\int\mathcal{D}\psi\mathcal{D}\bar{\psi}e^{-S_{0}}\left(-\delta
S+\frac{1}{2}\delta S^{2}+\ldots\right)\right]$ (82) $\displaystyle=\log
Z_{0}-\langle\delta S\rangle+\frac{1}{2}\langle\delta S^{2}\rangle+\ldots\,,$

where the expectation values are evaluated with the isotropic action. There
are two components of $\delta S$ that are $\order{\Phi^{2}}$ or lower, as can
be seen from Eq. (75). The first comes from the velocity anisotropy term
$\Phi\bar{\psi}\mu^{y}(\gamma^{x}\bar{D}_{x}-\gamma^{y}\bar{D}_{y})\psi$,
where the bar indicates that the covariant derivative is defined with the
isotropic monopole gauge configuration. The second component arises from
$\order{\Phi^{2}}$ corrections to the isotropic gauge configuration, which
appear in the term $\bar{\psi}\not{D}\psi$. Only the first of these
corrections gives $\order{\Phi^{2}}$ contributions to the free energy; the
corrections to the saddle-point gauge configuration $\delta A_{n}$ couple to
the conserved current $J^{\mu}=\bar{\psi}\gamma^{\mu}\psi+\order{\Phi}$, whose
expectation value vanishes. Hence, the free energy to $\order{\Phi^{2}}$ is
given by

$\displaystyle F$ $\displaystyle=\log
Z_{0}+2\Phi^{2}\Tr\int\differential[3]{r}\differential[3]{r^{\prime}}\langle\bar{\psi}(r)(\gamma^{x}(\bar{D}_{x}-\hat{x})-\gamma^{y}(\bar{D}_{y}-\hat{y}))\psi(r)$
(83)
$\displaystyle\times\bar{\psi}(r^{\prime})(\gamma^{x}(\bar{D}_{x}-\hat{x})-\gamma^{y}(\bar{D}_{y}-\hat{y}))\psi(r^{\prime})\rangle$

The additional factors of $\hat{x}\,,\hat{y}$ arise from the Weyl rescaling.

The calculation in Eq. (83) amounts to calculating the two-point function of
the fermion bilinear
$\bar{\psi}(r)(\gamma^{y}(\bar{D}_{x}-\hat{x})+\gamma^{x}(\bar{D}_{y}-\hat{y}))\psi(r)$.
To this end, we denote the matrix elements of the operator
$(\gamma^{y}(\bar{D}_{x}-\hat{x})+\gamma^{x}(\bar{D}_{y}-\hat{y}))$ in the
spinor harmonic basis

$B_{n,\ell
m\ell^{\prime}m^{\prime}}(\omega)=\int\differential{\Omega}\begin{pmatrix}T^{\dagger}_{n,\ell
m}(\theta,\phi)e^{i\omega\tau}\\\ S^{\dagger}_{n,\ell
m}(\theta,\phi)e^{i\omega\tau}\end{pmatrix}\left[\gamma^{y}(\bar{D}_{x}-\hat{x})+\gamma^{x}(\bar{D}_{y}-\hat{y})\right]\begin{pmatrix}T_{n,\ell
m}(\theta,\phi)e^{-i\omega\tau}&S_{n,\ell
m}(\theta,\phi)e^{-i\omega\tau}\end{pmatrix}$ (84)

These functions are exactly calculable in terms of Wigner $3$-$j$ symbols, and
are only non-zero for
$\absolutevalue{\ell-\ell^{\prime}}\,,\absolutevalue{m-m^{\prime}}\leq 2$. In
order to calculate these functions, we need the matrix elements

$\displaystyle\bra{Y_{q,\ell m}}\hat{x}\ket{Y_{q,\ell^{\prime}m^{\prime}}}$
(85) $\displaystyle\bra{Y_{q,\ell
m}}\hat{y}\ket{Y_{q,\ell^{\prime}m^{\prime}}}$ $\displaystyle\bra{Y_{q,\ell
m}}D^{\bot}_{x}\ket{Y_{q,\ell^{\prime}m^{\prime}}}$
$\displaystyle\bra{Y_{q,\ell
m}}D^{\bot}_{y}\ket{Y_{q,\ell^{\prime}m^{\prime}}}$

where $Y_{q,\ell m}$ is the scalar monopole harmonic [48] in a background
monopole of strength $2q\equiv n$ and $D^{\bot}_{i}$ is the angular component
of the covariant derivative $D_{i}$; the radial component is simply equal to
$\hat{i}\partialderivative{\tau}$. For this, we need the integral formula for
three monopole harmonics

$\displaystyle\int\differential{\hat{n}}Y_{q,\ell
m}Y_{q^{\prime},\ell^{\prime}m^{\prime}}Y_{q^{\prime\prime},\ell^{\prime\prime}m^{\prime\prime}}=(-1)^{\ell+\ell^{\prime}+\ell^{\prime\prime}}\sqrt{\frac{(2\ell+1)(2\ell^{\prime}+1)(2\ell^{\prime\prime}+1)}{4\pi}}\begin{pmatrix}\ell&\ell^{\prime}&\ell^{\prime\prime}\\\
q&q^{\prime}&q^{\prime\prime}\end{pmatrix}\begin{pmatrix}\ell&\ell^{\prime}&\ell^{\prime\prime}\\\
m&m^{\prime}&m^{\prime\prime}\end{pmatrix}$ (86)

The first two matrix elements can be easily computed with the identity

$\displaystyle\hat{x}$
$\displaystyle=-\sqrt{\frac{4\pi}{6}}\left(Y_{0,l,1}-Y_{0,1,-1}\right)$ (87)
$\displaystyle\hat{y}$
$\displaystyle=i\sqrt{\frac{4\pi}{6}}\left(Y_{0,l,1}+Y_{0,1,-1}\right)$
$\displaystyle\bra{Y_{q,\ell,m}}\hat{x}\ket{Y_{q,\ell^{\prime},m^{\prime}}}$
$\displaystyle=(-1)^{\ell+\ell^{\prime}+q+m}\sqrt{\frac{(2\ell+1)(2\ell^{\prime}+1)}{2}}\begin{pmatrix}\ell&1&\ell^{\prime}\\\
-q&0&q\end{pmatrix}\left[\begin{pmatrix}\ell&1&\ell^{\prime}\\\
-m&1&m^{\prime}\end{pmatrix}-\begin{pmatrix}\ell&1&\ell^{\prime}\\\
-m&-1&m^{\prime}\end{pmatrix}\right]$
$\displaystyle\bra{Y_{q,\ell,m}}\hat{y}\ket{Y_{q,\ell^{\prime},m^{\prime}}}$
$\displaystyle=-i(-1)^{\ell+\ell^{\prime}+q+m}\sqrt{\frac{(2\ell+1)(2\ell^{\prime}+1)}{2}}\begin{pmatrix}\ell&1&\ell^{\prime}\\\
-q&0&q\end{pmatrix}\left[\begin{pmatrix}\ell&1&\ell^{\prime}\\\
-m&1&m^{\prime}\end{pmatrix}+\begin{pmatrix}\ell&1&\ell^{\prime}\\\
-m&-1&m^{\prime}\end{pmatrix}\right]$

To calculate the last two matrix elements, we must utilize the raising and
lowering angular momenta operators

$\displaystyle
zL_{x}=z(L_{+}+L_{-})=D_{x}^{\bot}-\hat{y}\partialderivative{\phi}$ (88)

and similarly for $D_{y}^{\bot}$. This equation can be easily verified for the
angular momenta operators without a monopole background, and we verify
numerically that this formula correctly generalizes to non-zero $q$. This
leads to the formula

$\displaystyle\bra{Y_{q,\ell
m}}D^{\bot}_{x}\ket{Y_{q,\ell^{\prime}m^{\prime}}}=(-1)^{1+\ell+\ell^{\prime}+q+m}\sqrt{\frac{(2\ell+1)(2\ell^{\prime}+1)}{4}}\begin{pmatrix}\ell&1&\ell^{\prime}\\\
-q&0&q\end{pmatrix}$ (89)
$\displaystyle\times\left[\sqrt{(\ell^{\prime}-m^{\prime})(\ell^{\prime}+m^{\prime}+1)}\begin{pmatrix}\ell&1&\ell^{\prime}\\\
-m&0&m^{\prime}+1\end{pmatrix}-\sqrt{(\ell^{\prime}+m^{\prime})(\ell^{\prime}-m^{\prime}+1)}\begin{pmatrix}\ell&1&\ell^{\prime}\\\
-m&0&m^{\prime}-1\end{pmatrix}\right]$ $\displaystyle-$ $\displaystyle
im^{\prime}\bra{Y_{q,\ell,m}}\hat{y}\ket{Y_{q,\ell^{\prime},m^{\prime}}}$
$\displaystyle\bra{Y_{q,\ell
m}}D^{\bot}_{y}\ket{Y_{q,\ell^{\prime}m^{\prime}}}=-i(-1)^{1+\ell+\ell^{\prime}+q+m}\sqrt{\frac{(2\ell+1)(2\ell^{\prime}+1)}{4}}\begin{pmatrix}\ell&1&\ell^{\prime}\\\
-q&0&q\end{pmatrix}$
$\displaystyle\times\left[\sqrt{(\ell^{\prime}-m^{\prime})(\ell^{\prime}+m^{\prime}+1)}\begin{pmatrix}\ell&1&\ell^{\prime}\\\
-m&0&m^{\prime}+1\end{pmatrix}+\sqrt{(\ell^{\prime}+m^{\prime})(\ell^{\prime}-m^{\prime}+1)}\begin{pmatrix}\ell&1&\ell^{\prime}\\\
-m&0&m^{\prime}-1\end{pmatrix}\right]$ $\displaystyle-$ $\displaystyle
im^{\prime}\bra{Y_{q,\ell m}}\hat{y}\ket{Y_{q,\ell^{\prime}m^{\prime}}}$

From these matrix elements, the components of $B_{n,\ell
m\ell^{\prime}m^{\prime}}$ can be assembled by expressing the spinor monopole
harmonics in terms of the scalar harmonics.

Upon obtaining an expression for $B$, we have

$\displaystyle F=\log
Z_{0}-4\Phi^{2}\sum_{\ell,\ell^{\prime},m,m^{\prime}}\int\frac{\differential{\omega}}{2\pi}\Tr\left[B_{n,\ell
m\ell^{\prime}m^{\prime}}(\omega)\mathbf{N}_{n,\ell}(\omega+i\mathbf{M}_{n,\ell})B^{\dagger}_{n,\ell^{\prime}m^{\prime}\ell
m}(\omega)\mathbf{N}_{n,\ell^{\prime}}(\omega-i\mathbf{M}_{n,\ell^{\prime}})\right]$
(90)

The minus sign outside the summation relative to Eq. 83 arises from the
fermion loop. What remains is a suitable procedure for regularizing the
divergent expression in Eq. 90 \- as the functions $B_{n,\ell
m\ell^{\prime}m^{\prime}}$ are rather complicated summations of Wigner $3$-$j$
symbols, this is a non-trivial task and we leave this as an open question for
future study.

## References

  * Anderson [1973] P. Anderson, “Resonating valence bonds: A new kind of insulator?” Materials Research Bulletin 8, 153 (1973).
  * Zhou  _et al._ [2017] Y. Zhou, K. Kanoda, and T.-K. Ng, “Quantum spin liquid states,” Rev. Mod. Phys. 89, 025003 (2017).
  * Savary and Balents [2016] L. Savary and L. Balents, “Quantum spin liquids: a review,” Reports on Progress in Physics 80, 016502 (2016).
  * Ioffe and Larkin [1988] L. B. Ioffe and A. I. Larkin, “Effective action of a two-dimensional antiferromagnet,” Int. J. Mod. Phys. B 02, 203 (1988).
  * Gelfand  _et al._ [1989] M. P. Gelfand, R. R. P. Singh, and D. A. Huse, “Zero-temperature ordering in two-dimensional frustrated quantum Heisenberg antiferromagnets,” Phys. Rev. B 40, 10801 (1989).
  * Dagotto and Moreo [1989] E. Dagotto and A. Moreo, “Phase diagram of the frustrated spin-1/2 Heisenberg antiferromagnet in 2 dimensions,” Phys. Rev. Lett. 63, 2148 (1989).
  * Chandra  _et al._ [1990] P. Chandra, P. Coleman, and A. I. Larkin, “Ising transition in frustrated Heisenberg models,” Phys. Rev. Lett. 64, 88 (1990).
  * Read and Sachdev [1991] N. Read and S. Sachdev, “Large-$N$ expansion for frustrated quantum antiferromagnets,” Phys. Rev. Lett. 66, 1773 (1991).
  * Sachdev and Read [1991] S. Sachdev and N. Read, “Large $N$ Expansion for Frustrated and Doped Quantum Antiferromagnets,” Int. J. Mod. Phys. B 5, 219 (1991), arXiv:cond-mat/0402109 [cond-mat.str-el] .
  * Read and Sachdev [1989] N. Read and S. Sachdev, “Valence-bond and spin-Peierls ground states of low-dimensional quantum antiferromagnets,” Phys. Rev. Lett. 62, 1694 (1989).
  * Read and Sachdev [1990] N. Read and S. Sachdev, “Spin-Peierls, valence-bond solid, and Néel ground states of low-dimensional quantum antiferromagnets,” Phys. Rev. B 42, 4568 (1990).
  * Motrunich and Vishwanath [2004] O. I. Motrunich and A. Vishwanath, “Emergent photons and transitions in the O(3) sigma model with hedgehog suppression,” Phys. Rev. B 70, 075104 (2004), cond-mat/0311222 .
  * Senthil  _et al._ [2004a] T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P. A. Fisher, “Deconfined Quantum Critical Points,” Science 303, 1490 (2004a), arXiv:cond-mat/0311326 [cond-mat.str-el] .
  * Senthil  _et al._ [2004b] T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and M. P. A. Fisher, “Quantum criticality beyond the Landau-Ginzburg-Wilson paradigm,” Phys. Rev. B 70, 144407 (2004b), arXiv:cond-mat/0312617 [cond-mat.str-el] .
  * Liu  _et al._ [2021] W.-Y. Liu, J. Hasik, S.-S. Gong, D. Poilblanc, W.-Q. Chen, and Z.-C. Gu, “The emergence of gapless quantum spin liquid from deconfined quantum critical point,” (2021), arXiv:2110.11138 [cond-mat.str-el] .
  * Becca [2020] F. Becca, “Talk at KITP, University of California, Santa Barbara, November 5,” Correlated Systems with Multicomponent Local Hilbert Spaces (2020).
  * Wang and Sandvik [2018] L. Wang and A. W. Sandvik, “Critical Level Crossings and Gapless Spin Liquid in the Square-Lattice Spin-1 /2 $J_{1}$-$J_{2}$ Heisenberg Antiferromagnet,” Phys. Rev. Lett. 121, 107202 (2018), arXiv:1702.08197 [cond-mat.str-el] .
  * Ferrari and Becca [2020] F. Ferrari and F. Becca, “Gapless spin liquid and valence-bond solid in the $J_{1}$-$J_{2}$ Heisenberg model on the square lattice: Insights from singlet and triplet excitations,” Phys. Rev. B 102, 014417 (2020), arXiv:2005.12941 [cond-mat.str-el] .
  * Nomura and Imada [2020] Y. Nomura and M. Imada, “Dirac-type nodal spin liquid revealed by machine learning,” (2020), arXiv:2005.14142 [cond-mat.str-el] .
  * Liu  _et al._ [2020] W.-Y. Liu, S.-S. Gong, Y.-B. Li, D. Poilblanc, W.-Q. Chen, and Z.-C. Gu, “Gapless quantum spin liquid and global phase diagram of the spin-1/2 $J_{1}$-$J_{2}$ square antiferromagnetic Heisenberg model,” (2020), arXiv:2009.01821 [cond-mat.str-el] .
  * Capriotti  _et al._ [2001] L. Capriotti, F. Becca, A. Parola, and S. Sorella, “Resonating Valence Bond Wave Functions for Strongly Frustrated Spin Systems,” Phys. Rev. Lett. 87, 097201 (2001), arXiv:cond-mat/0107204 [cond-mat.str-el] .
  * Ivanov and Senthil [2002] D. A. Ivanov and T. Senthil, “Projected wave functions for fractionalized phases of quantum spin systems,” Phys. Rev. B 66, 115111 (2002), arXiv:cond-mat/0204043 [cond-mat.str-el] .
  * Hu  _et al._ [2013] W.-J. Hu, F. Becca, A. Parola, and S. Sorella, “Direct evidence for a gapless $\mathbb{Z}_{2}$ spin liquid by frustrating Néel antiferromagnetism,” Phys. Rev. B 88, 060402 (2013), arXiv:1304.2630 [cond-mat.str-el] .
  * Ferrari and Becca [2018] F. Ferrari and F. Becca, “Spectral signatures of fractionalization in the frustrated Heisenberg model on the square lattice,” Phys. Rev. B 98, 100405 (2018), arXiv:1805.09287 [cond-mat.str-el] .
  * Wen [1991] X. G. Wen, “Mean-field theory of spin-liquid states with finite energy gap and topological orders,” Phys. Rev. B 44, 2664 (1991).
  * Kitaev [2003] A. Kitaev, “Fault tolerant quantum computation by anyons,” Annals Phys. 303, 2 (2003), arXiv:quant-ph/9707021 .
  * Senthil and Fisher [2000] T. Senthil and M. P. A. Fisher, “$\mathbb{Z}_{2}$ gauge theory of electron fractionalization in strongly correlated systems,” Phys. Rev. B 62, 7850 (2000), cond-mat/9910224 .
  * Wen [2002] X.-G. Wen, “Quantum orders and symmetric spin liquids,” Phys. Rev. B 65, 165113 (2002), arXiv:cond-mat/0107071 [cond-mat.str-el] .
  * Senthil and Lee [2005] T. Senthil and P. A. Lee, “Cuprates as doped U (1) spin liquids,” Phys. Rev. B 71, 174515 (2005), arXiv:cond-mat/0406066 [cond-mat.str-el] .
  * Kitaev [2006] A. Kitaev, “Anyons in an exactly solved model and beyond,” Annals Phys. 321, 2 (2006), arXiv:cond-mat/0506438 .
  * Shackleton  _et al._ [2021] H. Shackleton, A. Thomson, and S. Sachdev, “Deconfined criticality and a gapless ${\mathbb{z}}_{2}$ spin liquid in the square-lattice antiferromagnet,” Phys. Rev. B 104, 045110 (2021).
  * Boyack  _et al._ [2018] R. Boyack, C.-H. Lin, N. Zerf, A. Rayyan, and J. Maciejko, “Transition between algebraic and ${\mathbb{z}}_{2}$ quantum spin liquids at large $n$,” Phys. Rev. B 98, 035137 (2018).
  * Dupuis  _et al._ [2021] E. Dupuis, R. Boyack, and W. Witczak-Krempa, “Anomalous dimensions of monopole operators at the transitions between Dirac and topological spin liquids,” (2021), arXiv:2108.05922 [cond-mat.str-el] .
  * Hermele  _et al._ [2005] M. Hermele, T. Senthil, and M. P. A. Fisher, “Algebraic spin liquid as the mother of many competing orders,” Phys. Rev. B 72, 104404 (2005), arXiv:cond-mat/0502215 [cond-mat.str-el] .
  * Franz  _et al._ [2002] M. Franz, Z. Tešanović, and O. Vafek, “${\mathrm{qed}}_{3}$ theory of pairing pseudogap in cuprates: From $d$-wave superconductor to antiferromagnet via an algebraic fermi liquid,” Phys. Rev. B 66, 054535 (2002).
  * Vafek  _et al._ [2002] O. Vafek, Z. Tešanović, and M. Franz, “Relativity restored: Dirac anisotropy in ${\mathrm{q}\mathrm{e}\mathrm{d}}_{3}$,” Phys. Rev. Lett. 89, 157003 (2002).
  * Borokhov  _et al._ [2002] V. Borokhov, A. Kapustin, and X. Wu, “Topological disorder operators in three-dimensional conformal field theory,” Journal of High Energy Physics 2002, 049 (2002).
  * Pufu [2014] S. S. Pufu, “Anomalous dimensions of monopole operators in three-dimensional quantum electrodynamics,” Phys. Rev. D 89, 065016 (2014).
  * Alicea [2008] J. Alicea, “Monopole quantum numbers in the staggered flux spin liquid,” Phys. Rev. B 78, 035126 (2008).
  * Song  _et al._ [2020] X.-Y. Song, Y.-C. He, A. Vishwanath, and C. Wang, “From spinon band topology to the symmetry quantum numbers of monopoles in Dirac spin liquids,” Phys. Rev. X 10, 011033 (2020), arXiv:1811.11182 [cond-mat.str-el] .
  * Song  _et al._ [2019] X.-Y. Song, C. Wang, A. Vishwanath, and Y.-C. He, “Unifying Description of Competing Orders in Two Dimensional Quantum Magnets,” Nature Commun. 10, 4254 (2019), arXiv:1811.11186 [cond-mat.str-el] .
  * Wang  _et al._ [2017] C. Wang, A. Nahum, M. A. Metlitski, C. Xu, and T. Senthil, “Deconfined quantum critical points: symmetries and dualities,” Phys. Rev. X 7, 031051 (2017), arXiv:1703.02426 [cond-mat.str-el] .
  * Huh and Sachdev [2008] Y. Huh and S. Sachdev, “Renormalization group theory of nematic ordering in $d$-wave superconductors,” Phys. Rev. B 78, 064512 (2008), arXiv:0806.0002 [cond-mat.str-el] .
  * Dupuis  _et al._ [2019] E. Dupuis, M. B. Paranjape, and W. Witczak-Krempa, “Transition from a dirac spin liquid to an antiferromagnet: Monopoles in a ${\mathrm{qed}}_{3}$-gross-neveu theory,” Phys. Rev. B 100, 094443 (2019).
  * Éric Dupuis and Witczak-Krempa [2021] Éric Dupuis and W. Witczak-Krempa, “Monopole hierarchy in transitions out of a dirac spin liquid,” Annals of Physics 435, 168496 (2021), special issue on Philip W. Anderson.
  * Thomson and Sachdev [2018] A. Thomson and S. Sachdev, “Fermionic Spinon Theory of Square Lattice Spin Liquids near the Néel State,” Phys. Rev. X 8, 011012 (2018), arXiv:1708.04626 [cond-mat.str-el] .
  * Metlitski  _et al._ [2008] M. A. Metlitski, M. Hermele, T. Senthil, and M. P. A. Fisher, “Monopoles in ${\text{cp}}^{N-1}$ model via the state-operator correspondence,” Phys. Rev. B 78, 214418 (2008).
  * Wu and Yang [1976] T. T. Wu and C. N. Yang, “Dirac monopole without strings: Monopole harmonics,” Nuclear Physics B 107, 365 (1976).