Anisotropic deconfined criticality in Dirac spin liquids

We analyze a Higgs transition from a 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\'eel 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\'eel 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 SU(2) spin liquid, which we have analyzed in prior work. By studying the deconfined critical theory separating the 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\'eel 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.

In a recent work with A. Thomson [31], we proposed a unified theory that contains multiple instabilities of the gapless Z 2 spin liquid Z2Azz13, 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 Z 2 spin liquid as a condensed phase of a parent SU (2) gauge theory coupled to Higgs bosons. This SU(2) theory, identified as the π-flux spin liquid, also has a proximate phase with 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 Z 2 spin liquid in the J 1 -J 2 model. In our previous work [31], we studied the transition between the SU(2) and Z 2 gauge theory using field-theoretic techniques; in this work, we complete our study by an analysis of the U(1) to 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 U(1) to 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 . The structure of our critical theory bears resemblence to prior studies of transitions between Dirac spin liquids and a gapped 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 [34][35][36] 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 Z 2 spin liquid of interest will lead to Lorentz symmetry breaking and a dynamical Towards a final phase diagram [1] L. Wang [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].

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 Z 2 spin liquid through a parent π-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 iα , α =↑ , ↓ at site i = (i x , i y ) of the square lattice using along with the constraints Introducing the Nambu spinor, and Pauli matrices τ which act on spinor indices, we can write a mean-field ansatz for our Hamil- where the hoppings u ij must be determined self-consistently. The additional degrees of freedom in this representation is reflected by an SU(2) g gauge symmetry, under which and a corresponding transformation for u ij . Including gauge fluctuations are necessary to enforce the constraint in Eq. (2.2), which in the Nambu spinor variables becomes ψ † i τ ψ i = 0. Different spin liquids may be described by different mean-field ansatzes u ij , which may also spontaneously break the relevant gauge fluctuations from SU(2) down to U(1) or Z 2 .
The particular spin liquid of relevance to the J 1 -J 2 model can be labeled as Z2Azz13 following Wen's classification [28]. This spin liquid, along with two relevant proximate spin liquids, can be described by the mean-field ansatz, The three relevant spin liquids are: • The π-flux phase with SU(2) gauge symmetry corresponds to χ = η = 0 and γ = 0.
• The Z 2 spin liquid Z2Azz13 is obtained from the staggered flux phase by turning on a non-zero γ.
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 π-flux phase have an emergent Lorentz invariance, whereas the staggered flux and Z 2 phase have anisotropic dispersion relations on a mean-field level (note that prior studies [34][35][36] 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 π-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 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 π-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 Z2Azz13 as well as these two instabilities can hence be obtained from QCD 3 with N = 2 fermion doublets and an SU(2) gauge group, and coupling it to Higgs fields whose condensation breaks the gauge group to either U(1) or 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. (2.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 ψ, an SU (2) gauge field A a µ , and three three-component adjoint Higgs fields Φ a 1,2,3 , a = x, y, z. The SU(2) gauge (σ) and valley (µ) Pauli matrices both rotate between the four fermion flavors.
We have absorbed the coefficients of the Yukawa couplings into the Higgs fields. The general form of the Higgs potential V (Φ) 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 I. The action of the SU(2) spin rotation symmetry requires a more careful analysis and is described below. All three Higgs fields transform trivially under this SU(2) symmetry. Our choice of representing the fermionic degrees of freedom in terms of Dirac fermions obfuscates the full SU(2) spin rotation symmetry, as rotations around the x or y axis involve charge conjugation. However, the U(1) subgroup corresponding to rotations around the z axis is simply given by a phase shift in ψ, ψ → e iθ ψ. The full SU(2) rotation is a π/2 rotation. We omit the action of SU(2) spin rotation symmetry as its action on the fermions ψ is non-trivial and described in the main text.
symmetry may be made explicit by writing the theory in terms of Majorana fermions [42], but this will not be necessary for our purposes.
On a mean-field level, this theory admits three phases, illustrated in Fig. 2. When all three Higgs fields are uncondensed, we recover the π-flux phase with SU(2) gauge group. When Φ 3 condenses, we obtain the U(1) staggered flux phase. The condensation of Φ 1 and Φ 2 yields the gapless spin liquid Z2Azz13. The masses of Φ 1 and Φ 2 are fixed to be equal by the microscopic square lattice symmetry, so both condense simultaneously. Furthermore, the symmetry-allowed cubic term abc Φ a 1 Φ b 2 Φ c 3 forces Φ 3 to condense along with Φ 1 , Φ 2 . Our conjectured trajectory of the J 1 -J 2 model as a function of J 2 J 1 is shown by the dotted blue line in Fig. 2, where transitions into either the staggered flux or π-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 π-flux phase; this deformation introduces first-order transitions near the tricritical point.
The transition from the π-flux phase to the Z 2 phase was studied in [31], and the focus of our work will be the U(1) to Z 2 transition. Before proceeding with our analysis, we briefly summarize the results of our study of the SU(2) to Z 2 transition. The primary order parameters for this theory are the masses of Φ 1,2 , so we neglect fluctuations of Φ 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 Φ 1,2 are divergent along a onedimensional 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 leadingorder 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 oneloop 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 −α exp β ln 2 r , with β = −6/π 2 for the VBS correlator and −12/π 2 for the Néel correlator, and α being some non-universal coefficient.
For the remainder of our work, we study the U(1) to 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 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 SU(4) symmetry and generate a non-zero velocity anisotropy, which further breaks the 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 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 O(N f ) effect on the scaling dimension of the monopole, in contrast to the Higgs screening, which is O(1).

III. LARGE N f EFFECTIVE ACTION
Our goal is to study the U(1) to Z 2 transition of Eq. (2.7). Both phases have the Higgs field Φ a 3 = 0, so we fix Φ a 3 = δ az Φ, Φ = 0. In this theory, the SU(2) gauge symmetry is broken to U(1), so we only consider A µ ≡ A 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Φ < 0, and the low energy behavior can be described in terms of a single complex Higgs This field transforms as a charge 2 Higgs field under the unbroken U(1) symmetry, as desired.
Other linear combinations of Φ 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 O(1/N f ) fixed point, which will lead to a first-order transition.
With these definitions, the Lagrangian that describes the U(1) → Z 2 transition is 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 Φ is a modification to the Lorentz-invariant fermion propagator 1/ / 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 Φ = 0 [34][35][36]. In our case, terms of this form are spontaneously generated at one-loop order by the critical Higgs field, and Φ 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 where the inverse propagators Γ, Π are generated by the one-loop fermion diagrams shown in  To calculate the effective propagators, we need the fermion propagator, which receives corrections to its Lorentz-invariant value of 1/ / p due to a non-zero Φ. This may be treated perturbatively in Φ, 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 Φ. We include further details of this calculation in Appendix C, and cite the results in the main text. Defining the variables the effective inverse Higgs propagator (obtained from the Φ-dependent fermion propagator) is Likewise, we need the general form of the effective gauge boson propagator. The presence of a non-zero Φ modifies the gauge coupling, and hence non-Lorentz-invariant corrections arise both from Φ-dependent modifications to the fermion propagator as well as O(Φ) vertices. We separate this calculation into three pieces. The first correction comes from using the O(Φ 0 ) vertices, but with the full fermion propagator. This one-loop term contributes The second correction comes from using one O(Φ) vertex, which gives the contribution There is also an extra factor of 2 in Π (2) xx,yy due to the two possible vertex orderings. Finally, the third correction comes from using two O(Φ) vertices, We verify that the combined inverse propagator Π µν (k) = i=1,2,3 Π (i) µν (k) annihilates the vector (k 0 , k x , k y ), as required by gauge invariance. Note that Π µν requires a gauge fixing term in order to be invertable. Followin Ref. [34], we add the following non-local gauge fixing term to the All gauge-invariant observables have been checked to ensure they are independent of the choice of ξ.

IV. RENORMALIZATION GROUP ANALYSIS
We perform a renormalization group (RG) analysis of the O(1/N f ) effective theory. We are interested in studying the behavior of this theory under the rescaling We also define a rescaling of the fermion fields 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.

A. Fermion self-energy
We first evaluate the O(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 Λ.
The logarithmic derivative of the fermion self-energy with respect to this cutoff takes the general 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 where F and G are homogeneous functions of the three-momenta, and K(y) serves as a UV cutoff with the property that K(0) = 1 and 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 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 wherep ≡ (cos θ, sin θ sin φ, sin θ cos φ). The resulting integrals in Eq. (4.5) are fully convergent and may be evaluated numerically, from which we can extract the coefficients C 0,1,2 .

B. Fixed points
The RG equations for the velocity anisotropy Φ are In the absence of the Higgs field, Φ has a stable fixed point at Φ = 0. The gauge field contribution to this equation has been calculated to leading order in Φ [34], and we verify agreement with this result. The evaluation of Eq. (4.6) is plotted in Fig. 4. A stable fixed point is found at Φ c ≈ 0.45765.
At this point, the dynamical critical exponent z is given by Recall that when we derived this critical U(1) → Z 2 theory as a component of a parent SU (2) theory, we made a gauge choice such that wΦ < 0, where w is the coefficient of the symmetry- When Φ a 3 = Φδ az , we can diagonalize this term to yield a mass, wΦ(H * H − M * M), where H is the combination of Φ x,y 1,2 given in Eq. 3.1 and M is a charge-2 Higgs field of a similar form but with x ↔ y. If we assume w > 0, Φ < 0, then H will become massless first, but the fixed-point value of Φ gives a negative mass to M, leading to a first-order transition driven by the condensation of M. As a consequence, we must fix our parent SU(2) theory to have w < 0 in order to yield a continuous transition. If we had made a gauge choice such that wΦ > 0, then our theory would have been driven by the condensation of M rather than H; this still leads to the gapless Z 2 spin liquid Z2Azz13, and all gauge-invariant observables at the critical point remain the same, although the sign of Φ c changes.

C. 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ψµ z,x ψ. As mentioned previously, our particular representation obfuscates the full SU (2) action of spin rotation symmetry; however, the U(1) subgroup generated by rotations around the z-axis is given by the global U(1) symmetry ψ → e iθ ψ (recall that this is not the U(1) gauge symmetry, which acts as e iθσ z ). As a consequence, we focus on the z-component of the Néel order parameter, which is given byψµ y ψ. 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 µ 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.

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 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 µ x,z (VBS order parameter) or µ y (Néel order parameter). Once again, we can take the logarithmic derivative and remove the explicit cutoff dependence, leading to the equation 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 Φ c gives (4.10) The Néel and VBS correlators in momentum space have the scaling form (4.11) Making a Fourier transform to real space, the equal-time Néel and VBS correlators have the power law decay 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 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.

D. 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 SU(4) flavor symmetry. We refer to Ref. [34] for a more extensive study of the intertwining physical order parameters of the SU(4) theory and which relations hold in the presence of the velocity anisotropy -for our purposes, we note that the emergent SO(5) ⊂ SU(4) symmetry that relates the Néel and VBS order parameters is broken down to the microscopic SO(3) × 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 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 Φ. The VBS correlator is more difficult to study non-perturbatively in the velocity anisotropy, and we will later compute corrections to leading order in Φ.
The two-point function in momentum space is given by, with Q(p) the fermion propagator, p 0 (p 0 + k 0 ) + ap x (p x + k x,a ) + ap y (p y + k y,a ) p 2 (p + k ± ) 2 As before, we define k ± = (k 0 , k x ± Φk y , k y ± Φk x ). The Fourier transform d 3 k (2π) 3 e ik·r |k ± | (4.14) can be performed by a change of variables to give where we change to spherical coordinates, t = r cos θ, x = r sin θ sin φ, y = r sin θ cos φ, and Therefore, Fig. 6. We note the enhanced correlations of the Néel order parameter along the diagonals, which holds true for generic values of Φ.
An analogous calculation of the VBS order parameter is less analytically tractable, as the oneloop integral cannot be made isotropic by a coordinate transformation. As such, we resort to a perturbative study of the velocity anisotropy. This gives 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 mostrelevant 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 µ 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 Z 2 spin liquid Z2Azz13 from the 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 Z 2 spin liquid may emerge as a Higgsed phase from either a U(1) or 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 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 SU(2) → Z 2 and U(1) → 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 SU(2) → 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 = 1), or lack of a full 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  The phase diagram of the theory we have defined contains three phases on a mean-field levelthe SU(2) spin liquid of the π-flux phase, the U(1) staggered flux phase, and the gapless Z 2 spin liquid whose projective symmetry group labels it as Z2Azz13 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 Z 2 spin liquids are accessible in this general framework by a modification of the Higgs couplings [46].
Although there are many Z 2 spin liquids accessible starting from the π-flux phase and we will not attempt to study all of them, we will identify the continuum theories associated with the eight Z 2 spin liquids proximate to both the π-flux and staggered flux phase.

Symmetry fractionalization in the continuum staggered flux theory
The 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 SU(2) gauge transformation. This concept holds true in our continuum theory, and we will use this to identify possible 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 π-flux phase to the staggered flux phase, which consists of four Dirac fermions and a three-component adjoint Higgs field Φ 3 , both minimally coupled to an SU(2) gauge field, and an additional Yukawa coupling When Φ 3 is uncondensed, it can be integrated out, and the Yukawa couplings generate terms irrelevant at long distances. When Φ 3 is condensed -for concreteness, Φ a 3 = Φδ az -we replace Φ a 3 by its expectation value, and the fermion bilinear ψσ z µ y (γ x i∂ y + γ y i∂ x )ψ naively breaks translation and rotation symmetry according the symmetry transformations in Table I. However, the fermion bilinear is invariant under a combination of microscopic symmetries and SU(2) gauge transformations: U G is the action of the microscopic symmetry and V G is a gauge transformation, given by where g(φ) ≡ e iφσz reflects the residual 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 π/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 I. The ability to condense additional Higgs fields to yield a Z 2 spin liquid is constrained by the requirement that there exists a choice of phases φ 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 φ 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.

Translationally-invariant spin liquid ansatzes
There are four proximate 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 meanfield ansatz. These spin liquids are Z2Azz13, Z2A001n, Z2Azz1n, and Z2A0013. These spin liquids are classified based off of their symmetry fractionalization; in other words, how symmetry operations like T −1 y T x T y T −1 x are not directly equivalent to the identity, but equivalent up to a gauge transformation. These values are shown in Table II, along with the values corresponding to the continuum U(1) staggered flux phase, which can be read off from the transformations A3. Using this table to identify the phases φ i we find the continuum projective symmetry group for Z2Azz13, Z2Azz1n, and Z2A001n, From this, we can identify the Yukawa couplings consistent with these transformations, given in Table III. In order to realize these transitions, we require two Higgs fields, Φ 1 and Φ 2 which transform into each other under square lattice rotations and condense to acquire a non-zero Φ x 1 and Φ y 2 . 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 Z2Azz13 transition.

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 II aside from a change in sign in rows 1, 7, and 8. This leads to the continuum projective symmetry group for Z2Bzz13, Z2Bzz1n, and Z2B001n, Describing these phases as a Higgsed phase of the π-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 Φ 1 , which acquires a non-zero expectation value Φ x 1 in the Z 2 phase. This prevents a direct transition from the SU(2) π-flux state to the Z 2 spin liquid, as the masses of two Higgs fields Φ 1 and Φ 3 aren't constrained to be equal by the microscopic symmetries.
Appendix B: Emergent subsystem symmetries in the SU(2) → Z 2 transition In this appendix, we provide a more detailed discussion of the emergent subsystem symmetry in the SU(2) → Z 2 transition present in our theory. This transition is driven by the simultaneous condensation of the two Higgs fields Φ 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 V (Φ) contains various symmetry-allowed potential terms for Φ 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 Φ 1 Yukawa coupling in invariant under the transformation and likewise for Φ 2 , with U (x) the SO(3) rotation corresponding to the adjoint SU(2) action of e ifa(x)σ a , i.e.
and likewise forŨ (y). The simplest way to see the existence of these symmetries is to consider Φ 1 and Φ 2 as the x and y components, respectively, of fictitious gauge fields corresponding to the SU(2) symmetries ψ → e iσ a µ z,x ψ. 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 Φ 1 (Φ 2 ) Yukawa coupling and the fermion quadratic term are invariant under x-dependent (y-dependent) 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 µ as well as the Higgs potential V (Φ); 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 Φ a 1,2 ∂ 2 Φ a 1,2 , and hence our analysis suggests that these emergent subsystem symmetries control the behavior of physical observables at criticality. The contributions to the Néel and VBS order parameters have logarithm squared divergences, which we conjecture should lead to correlations decaying as 1 r α e −β ln 2 (r) , with β = − 6 π 2 N f for the VBS correlator and − 12 π 2 N f for the Néel correlator, and α being some non-universal coefficient. Interestingly, we note that the logarithm squared divergences for correlations of the scalar spin chirality exactly cancel at O N −1 f . We expect that this is related to the invariance of the scalar spin chirality order parameter under the subsystem symmetries given by Eqs. (B2) and (B3).
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 SU(2) space given by the σ a operators, and as such, can rotate gauge-invariant operators such asψµ i ψ 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. (3.2), the bare fermion action is given bȳ It is useful to work in an eigenbasis of σ z µ y . Inverting this, we get a diagonal 4 × 4 matrix in gauge/valley space, with elements As in the main text, we define the variables The inverse Higgs propagator is generated by the one-loop diagram given by Fig. 3 of the main text, which translates to the integral 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 2 a (k + p) 2 a (C4) To simplify this integral, we perform a change of integration variables to (p 0 , p x,± , p y,± ). Performing this change of variables gives a factor of (1 − Φ 2 ) −1 in the integral.
To evaluate this integral, we calculate the two general forms of integrals relevant to Eq. (C4).
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. (3.5) for general p i , k i , we get the Higgs propagator 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 Σ(k) are UV divergent, and hence require a UV cutoff Λ. The behavior of the self-energy upon integrating out high-energy modes is dictated by the logarithmic derivative with respect to the cutoff, Λ d dΛ Σ(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 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 K(y) serves as a UV cutoff with the property that K(0) = 1 and K(y) falls off rapidly for large y, i.e., K(y) = e −y . Since we are interested in the behavior at small momenta, we expand around k = 0, We then take the logarithmic derivative, We now convert to spherical coordinates, p = yΛ(cos θ, sin θ sin φ, sin θ cos φ), and use the homogeneity property of F and G to pull out factors of (yΛ) −1 .
The integral over y can be done explicitly via integration by parts, which causes the dependence on the cutoff function K to drop out. This leads to the expression cited in the main text wherep ≡ (cos θ, sin θ sin φ, sin θ cos φ).
Explicitly, we take, for the Higgs contribution to the self-energy, defining Q(p) as the fermion propagator, and perform a change of variables to shift the anisotropy to the spatial coordinates which yields the real space correlator given in the main text.
To compute the Fourier transform of the VBS correlator, given perturbatively by we define the function The Fourier transform of this function can be calculated by a similar change of variables, The various terms in the O(Φ 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.
leading order in Φ, the gauge field configuration corresponding to the isotropic Dirac monopole will be sufficient.
We start with the action for QED 3 with the allowed velocity anisotropy term, omitting the Higgs fields as they will not play any role in the calculation 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 µ ≡ ∂ µ − iA µ , and the microscopic SU(2) spin rotation symmetry is implemented explicitly by the σ 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 bȳ The non-zero anisotropy will affect these saddle-point solutions. The leading order corrections to In order to calculate the scaling dimension of this monopole, we set r = e τ and perform a Weyl rescaling g µν → e −2τ g µν 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 ψ in terms of these harmonics, where T n, m , S n, m are eigenvalues of the orbital angular momentum operator L 2 in the presence of a strength n monopole, with orbital angular momentum and total angular momentum + 1/2 for T n, m and − 1/2 for S n, m . Explicit expressions for the spinor harmonics T n, m and S n, m may be found in [37]. The variables Ψ m T , Ψ m S are anti-commuting coefficients. Expanded in this form, the isotropic action with Φ = 0 can be written as where Y q, m is the scalar monopole harmonic [48] in a background monopole of strength 2q ≡ n and D ⊥ i is the angular component of the covariant derivative D i ; the radial component is simply equal toî ∂ ∂τ . For this, we need the integral formula for three monopole harmonics dn Y q, m Y q , m Y q , m = (−1) + + (2 + 1)(2 + 1)(2 + 1) 4π   q q q To calculate the last two matrix elements, we must utilize the raising and lowering angular momenta and similarly for D ⊥ y . 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 The minus sign outside the summation relative to Eq. F9 arises from the fermion loop. What remains is a suitable procedure for regularizing the divergent expression in Eq. F16 -as the functions B n, m m 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.