Anomalous dimensions of monopole operators at the transitions between Dirac and topological spin liquids

Monopole operators are studied in a large family of quantum critical points between Dirac and topological quantum spin liquids (QSLs): chiral and Z$_{2}$ QSLs. These quantum phase transitions are described by conformal field theories (CFTs): quantum electrodynamics in 2+1 dimensions with 2N flavors of two-component massless Dirac fermions and a four-fermion interaction. For the transition to a chiral spin liquid, it is the Gross-Neveu interaction (QED$_{3}$-GN), while for the transitions to Z$_{2}$ QSLs it is a superconducting pairing term with general spin/valley structure (generalized QED$_{3}$-Z$_{2}$GN). Using the state-operator correspondence, we obtain monopole scaling dimensions to sub-leading order in 1/N. For monopoles with a minimal topological charge q=1/2, the scaling dimension is 2N*0.26510 at leading-order, with the quantum correction being 0.118911(7) for the chiral spin liquid, and 0.102846(9) for the simplest Z$_{2}$ case (the expression is also given for a general pairing term). Although these two anomalous dimensions are nearly equal, the underlying quantum fluctuations possess distinct origins. The analogous result in QED$_{3}$ is also obtained and we find a sub-leading contribution of -0.038138(5), which differs slightly from the value first obtained in the literature. The scaling dimension of a QED$_{3}$-GN monopole with minimal charge is very close to the scaling dimensions of other operators predicted to be equal by a conjectured duality between QED$_{3}$-GN with 2N=2 flavors and the CP$^{1}$ model. Additionally, non-minimally charged monopoles on both sides of the duality have similar scaling dimensions. By studying the large-q asymptotics of the scaling dimensions in QED$_{3}$, QED$_{3}$-GN, and QED$_{3}$-Z$_{2}$GN we verify that the constant O(q$^{0}$) coefficient precisely matches the universal non-perturbative prediction for CFTs with a global U(1) symmetry.

Gauge theories play an important role in modern condensed matter physics, in part due to their ability to provide a low-energy description of many quantum phases of matter. Gauge fields emerge as collective excitations that capture the highly entangled nature of certain strongly correlated systems. This is notably apparent in the case of frustrated two-dimensional magnets hosting quantum spin liquids and deconfined quantum critical points (dQCPs).
In these lattice systems, the emergent gauge field is compact, and, as a result, has topological excitations created by topological disorder operators. For a U(1) gauge field, these objects are called monopole or instanton operators and they play an essential role in many physical systems. Crucially, monopole proliferation confines the gauge field. This is the case in the pure U(1) gauge theory [1; 2]. In the presence of massless matter, however, monopoles are screened and confinement can be avoided if enough flavors of massless matter are present.
In particular, we will be first considering a transition from a U(1) Dirac spin liquid (DSL), which is described by QED 3 with 2N flavors of massless two-component Dirac fermions. Realizations of the U(1) DSL were formulated for the Kagome Heisenberg spin-1/2 magnet [3][4][5][6][7][8] and the J 1 − J 2 spin-1/2 model on the triangular lattice [9][10][11] with 2N = 4 flavors. Dirac spin-orbital liquid with effective spin j = 3/2 and 2N = 8 flavors have also been formulated for quantum magnets on honeycomb [12; 13] and triangular [14] lattices. For a large number of fermion flavors 2N , it has been shown through a 1/N expansion that monopole operators are irrelevant [15], and thus the U(1) DSL is stable in this limit. Taking into account next-to-leading order corrections [16], the critical number of fermion flavors was estimated to be 2N c = 12, beyond which minimally charged monopoles become irrelevant. This result was confirmed by Monte Carlo computations [17] and is consistent with conformal bootstrap bounds [18; 19]. The addition of disorder renders the model more unstable [20].
Monopole operators also serve as order parameters in neighboring phases. For instance, in the CP 1 model, which describes the transition between an antiferromagnetic (AFM) phase and a valence bond solid (VBS), there are monopoles with lattice quantum numbers and their condensation results in VBS order [21][22][23][24]. It is in this model where the scaling dimension of monopole operators were first obtained [25]. Monopoles are also crucial in the U(1) DSL, a parent state for many spin liquids. In this fermionic theory, monopoles can carry different quantum numbers due to the existence of fermion zero modes [26] which may dress monopoles in various ways. Monopoles describe various VBS and AFM orders, depending on which lattice the U(1) DSL is formulated [27].
By tuning a flavor-dependent Gross-Neveu (GN) interaction, a fermion mass is generated and monopoles with specific quantum numbers condense [28]. In particular, the U(1) DSL on the Kagome lattice orders to an antiferromagnetic 120 • coplanar order as monopoles dressed with a magnetic spin polarization condense [29; 30]. This confinement-deconfinement transition is described by the QED 3 -chiral Heisenberg GN model (QED 3 -cHGN), and the scaling dimensions of monopole operators at the quantum critical point (QCP) were obtained in Ref. [31]. In this case, the activation of the cHGN interaction breaks the flavor symmetry, resulting in a hierarchy among monopoles where the scaling dimension depends on the total magnetic spin of a monopole [32; 33].
The quantum criticality of the dQCP and the U(1) DSL with 2N = 4 fermions were recently given a precise relation. It was shown that they can be formulated as so-called Stiefel liquids, which are related to non-linear sigma models in 2 + 1 dimensions with target manifolds SO (n) / SO (4), where n = 5 and n = 6 for the dQCP and the U(1) DSL respectively [34]. Higher values are conjectured to realize non-lagrangian critical systems, for instance realizing a phase between a non-coplanar magnet and a VBS order when n = 7.
In this work, we will focus on transitions from the U(1) DSL to two topological quantum spin liquids (QSL): a chiral spin liquid (CSL) as mentioned above, and a general type of Z 2 QSL. The first transition is described by QED 3 -GN [35; 36], where the CSL results from the condensation of a symmetric fermion mass induced by the GN interaction. This transition can be realized for the Kagome [37][38][39][40] and triangular [41][42][43] Heisenberg magnets with 2N = 4 Dirac cones. The QCP for the noncompact QED 3 -GN has been studied in Refs. [36; 44-48]. Even though a mass gap is condensed in the CSL, there is no confinement-deconfinement transition taking place in the compact theory. Despite removing the screening effect of gapless modes, the symmetric condensed mass induces a Chern-Simons term in the infrared limit which gaps the monopoles and prevents their proliferation. The spinons in the gapped CSL thus remain deconfined. The CSL is a topologically ordered state that breaks timereversal symmetry, and has robust chiral edge modes.
In contrast, the non-chiral Z 2 QSL is obtained when the fermionic spinons undergo a pairing instability to a gapped s-wave superconducting state. The U(1) gauge field is gapped through the Higgs mechanism, and gives place to a discrete Z 2 gauge field. In this case, fractionalization remains intact. For the simplest case where the pairing interaction is diagonal in flavor space, the corresponding quantum phase transition was studied in Refs. [49; 50]. Earlier studies [30; 51] qualitatively described how a Z 2 QSL can be obtained from a Dirac QSL through a superconducting transition for the fermions, albeit without a fluctuating scalar field (Cooper pair field). In addition, Ref. [52] studied a similar model in the context of superconducting criticality in topological insulators. Interestingly, it turns out that, at leading-order in 1/N , the monopoles have the same scaling dimension at both QCPs as in the U(1) DSL [31; 50]. In this work, we obtain the next-to-leading order correction to monopole scaling dimensions at those QCPs. Futhermore, we also consider a more general class of Z 2 QSLs where the pairing interaction is not the same for all spin and valley degrees of freedom. We compute the anomalous dimension for the general Z 2 QSLs, and we also determine their bandstructure and Chern number inside the spin liquid phase.
This study is also motivated by the duality between QED 3 -GN with 2N = 2 fermion flavors and CP N −1 with N = 2 complex boson flavors conjectured in Ref. [53] and further studied in Ref. [44]. This duality can be checked by comparing the scaling of monopole operators with various scaling dimensions that are predicted to be equal according to this duality. The good agreement obtained in the LO result [31] is further improved by the scaling dimension correction we obtain here for the QED 3 -GN monopoles.
The paper is organized as follows. In the next section, we present the QED 3 -GN model and show how the stateoperator correspondence is used to obtain monopole scaling dimensions. In Sec. III, the leading-order computation presented in Ref. [31] is reviewed. In Sec. IV, 1/N corrections to monopole scaling dimensions are computed. We also verify that the scaling dimensions satisfy a conjectured convexity property. In Sec. V, we compare our results with the large-charge expansion obtained in Ref. [54] for CFTs. In Sec. VI, we study the QED 3 -GN ⇔ CP 1 duality [53]. In Sec. VII, we study monopole scaling dimensions at the transition to a Z 2 QSL, and obtain distinct values compared to the CSL. In Sec. VIII, we briefly discuss other phase transitions that could be studied with this formalism, including the QED 3 -U(N ) × U(N ) GN, QED 3 -chiral XY GN, and QED 3 -cHGN QCPs. We conclude with a discussion of our results and an outlook. In Appendix A, we review the phase transition from the U(1) DSL to the CSL in the non-compact model. In Apps. B and C, we give more details regarding how the kernels appearing in Sec. IV are obtained and simplified with gauge invariance. The expansion of these kernels in terms of harmonics is detailed in Apps. D and E. We give detailed simplifications of the kernels used for the case of minimally charge monopoles in App. F. In App. G, some remainder coefficients used to analytically approximate sums over angular momenta are shown. In App. H, we show how some contributions of fermion zero modes neglected in the main text vanish. In App. I, the fitting procedure used to alleviate finite-size effects when computing monopole anomalous dimensions are described. In App. J, we list monopole anomalous dimensions in QED 3 , QED 3 -GN, and QED 3 -Z 2 GN for topological charges up to q = 13.

II. MONOPOLES AT TRANSITION BETWEEN DIRAC & CHIRAL SPIN LIQUIDS
The action of the QED 3 -GN model in euclidean flat spacetime is given by (2.1) where Ψ is a 2N flavor spinor Ψ = (ψ 1 , ψ 2 , . . . ψ 2N ) with each flavor ψ i being a two-component Dirac fermion. For certain quantum magnets, where fermions emerge as fractionalized quasiparticles, the 2N flavors are related to two magnetic spin polarizations s =↑, ↓ and N valley nodes per spin, v = 1, . . . , N . Typical quantum magnets have N = 2 or 4 nodes, but here we keep N general and use it as an expansion parameter. The adjoint spinor is given by Ψ = Ψ † γ 0 , where the gamma matrices are defined in terms of the Pauli matrices by γ x,y = σ x,y , and γ 0 = σ z . The fermions are coupled to a compact U(1) gauge field A µ , and have a GN self-interaction with coupling strength h. The ellipsis denotes an irrelevant Maxwell term and the contribution of monopole operators M q (x) that we discuss further in what follows.
In 2 + 1 dimensions, U(1) gauge theories have an extra global U top (1) symmetry associated with the following conserved current where "top" stands for topological. The operators charged under U top (1) are called topological disorder operators or instantons. In this 2 + 1 dimensional context, we refer to them as monopole operators. These operators create topological configurations of the gauge field A q µ with a quantized flux dn µ µνρ ∂ ν A q ρ = 4πq, where the topological charge is a half-integer q ∈ Z/2 as a result of the Dirac quantization condition [55]. These kinds of configurations are allowed in the compact formulation of the U(1) gauge group, which gives the correct description for emergent gauge theories in a condensed matter context. The monopole operators themselves can be defined by the action of the topological current on them: where the ellipsis denotes less singular terms in the operator-product expansion (OPE) [15]. The resulting factor in front of the monopole operator corresponds to the magnetic field of a charge-q Dirac magnetic monopole. The model in Eq. (2.1) describes a transition from a DSL to a CSL. For a sufficiently strong coupling, a chiral order develops due to the condensation of a fermion bilinear: ΨΨ = 0. This may be studied by introducing an auxiliary pseudo-scalar boson φ. The effective action at the quantum critical point (QCP), denoted by S c eff , is where φ is an auxiliary boson decoupling the GN interaction. More details are shown in App. A.
In the compact version of QED 3 -GN, monopole operators are also present at the QCP. The main goal of this paper is to compute their scaling dimension ∆ Mq , which controls the scaling of the monopole two-point correlation function: Since the QED 3 -GN model at the QCP is a conformal field theory (CFT), the state-operator correspondence can be used to obtain these scaling dimensions [15]. This correspondence relies on a radial quantization of the CFT and a conformal transformation mapping the dilatation operatorD on R 3 to a HamiltonianĤ on S 2 × R. Denoting the usual radius on R 3 as r = e τ , 1 the related Weyl transformation of the spacetime is written as The scaling dimension of an operator O(x) then corresponds to the energy of some stateĤ |O = ∆ O |O on this compactified spacetime. Specifically, the charge-q operator with the smallest scaling dimension corresponds to the ground state of the CFT on the compactified spacetime S 2 × R, where the sphere S 2 is pierced by 4πq flux. To implement this flux, an external gauge field is coupled to the fermions or A q φ = (1 − cos θ) / sin θ in component notation. The smallest scaling dimension of monopole operators in topological sector q is then given by 2 where F q is the free energy and Z [A q ] is the partition function formulated on S 2 × S 1 β , i.e., the previous spacetime but now with the "time" direction compactified to a "thermal" circle S 1 β with radius β. This formulation allows us to introduce the holonomy of the gauge field along this circle, written as The holonomy couples to the fermion number operator d 2 r g(r)Ψ † Ψ =N fermions and acts as a chemical potential [19; 31]. The saddle-point equation of this holonomy constrains the fermion number to vanish = N fermions , (2.10) 1 We work in natural units where the two-sphere radius is R = 1. 2 More formally, we could write ∆q = lim β→∞ (Fq − F 0 ) as explained in Ref. [56]. However, it turns out that lim β→∞ F 0 = 0.
where "s.p." stands for saddle-point. The holonomy thus serves as a Lagrange multiplier that ensures that a state with N fermions = 0 is selected to correctly represent a gauge-invariant monopole operator [31].
The scaling dimension will be obtained using a large-N expansion. We first note that the partition function can be written as a path integral: (2.11) The effective action is now given by where / D A+A q is the gauge-covariant derivative on a curved spacetime including the external gauge field A q µ sourcing the 4πq flux: The gamma matrices γ b still correspond to the Pauli matrices, as the spacetime index is normalized with a tetrad e µ b which encapsulates the information about the metric g µν e µ b e ν c = δ bc . The path integral defining the partition function can be expanded around the saddle-point values of the auxiliary and gauge bosons: 14) which are defined by the following saddle-point conditions δF q δφ φ= φ ,Aµ= Aµ = δF q δA µ φ= φ ,Aµ= Aµ = 0. (2.15) Taking the fluctuations to scale as 1/ √ 2N , the saddlepoint expansion of the partition function is then where S (2) eff is the second variation of the action. Integrating over the quadratic fluctuations, this gives us the 1/N expansion of the free energy: . (2.18) Using the relation in Eq. (2.8), which follows from the state-operator correspondence, these first two terms of the free energy give the scaling dimension at next-toleading order in 1/N . 3 3 The expansion is in terms of the total number of fermion flavors, Since the fermionic mass condensed in the ordered phase is flavor-symmetric ΨΨ , the global flavor symmetry remains unbroken and monopole operators are organized as representations of SU(2N ). Just as for the various fermion bilinears and monopole correlation functions in the U(1) DSL [29; 57; 58], monopole correlation functions at the QCP between U(1) DSL and CSL related by this SU(2N ) symmetry share the same scaling dimension. 4 Depending on the lattice, various magnetic and VBS correlation functions will be described by minimally charged monopole operators [27], but they all share the same scaling dimension 2N × 0.26510 + 0.118911(7) + O N −1 , where the leading-order was found in Ref. [31] and the next-to-leading order is one of the main results of this work shown in Eq. (4.61). For typical quantum magnets, we have 2N = 4 fermion flavors. The way that monopole scaling dimensions control observable correlation functions could also be compared at this QCP and deep in the U(1) DSL phase. In this latter case, scaling dimensions are those of monopoles in QED 3 .

III. REVIEW OF N = ∞ THEORY
First, we review the computation of monopole scaling dimensions in QED 3 -GN at leading-order in 1/N [31]. At this order, the free energy is given by the effective action in (2.12) at its saddle-point corresponding to a global minimum where the trace over the 2N flavors has been taken and has canceled a prefactor of (2N ) −1 . The expectation value of the pseudo-scalar field is taken to be homogeneous. The gauge field is also constant at the saddlepoint, with a possible non-vanishing holonomy α on the thermal circle described in (2.9). The determinant operator is diagonalized by introducing monopole harmonics Y q, ,m (n), which are a generalization of spherical harmonics for a space with a charge at the center [59; 60]. For a fixed charge q, these functions form a complete basis. One important difference with these harmonics is that their angular momentum is now bounded below by this charge q. Using these functions to build appropriate eigenspinors, the eigenvalues of this determinant opera-4 This degeneracy we described is among a flavor symmetry multiplet composed of monopoles with the smallest scaling dimension, which we focus on. We emphasize that monopoles with larger scaling dimensions can be built by dressing fermion modes with higher energies. For instance, a splitting of monopoles was obtained for QED 3 q = 1 monopoles as their scaling dimensions increases with Lorentz spin [19]. Notably, for a monopole with Lorentz spin of order √ N , there is a O(N 0 ) additional positive correction.
tor on S 2 × S 1 β in Eq. (3.1) are shown to be [15; 31] −i × ω n − iα + iε q = q, where, for simplicity, we suppose q > 0 throughout. Here, ω n = 2πβ −1 (n + 1/2), for n ∈ Z, are the fermionic Matsubara frequencies, and ε are the energies of the modes for the quantized theory on S 2 × R: More details on the diagonalization are presented in App. (D 2). Note that the energies are dimensionless, as we work in units where the radius of the sphere is 1.
Each mode has the usual degeneracy that comes from the azimuthal symmetry, d = 2 . The free energy at leading-order then becomes The saddle-point equation for the holonomy, given in Eq. (2.15), yields the condition 2d sinh (βα) cosh (βε ) + cosh (βα) = 0, (3.5) which is solved for α = φ in the β → ∞ limit. With this result, the second gap equation at leading-order in β is given by whose only solution is φ = 0. Therefore, the saddlepoint values of both fields vanish. Inserting this result in Eq. (3.4), the monopole scaling dimension at leading-order in 1/N is obtained from Eq. (2.8) 5 where the energies at the saddle-point are defined as (3.8) 5 More explicitly, the leading-order free energy is lim β→∞ F Using zeta regularization, it follows that the q = 0 case vanishes: lim β→∞ F  This is simply the leading-order scaling dimension of QED 3 [15] (which must still be regularized). For example, the scaling dimension of the monopole with minimal charge is ∆ q=1/2 = 2N × 0.265 + O N 0 . Here, a supplementary GN interaction is considered, but it does not come into play at this level of the expansion since φ = 0. Thus, monopoles in QED 3 and QED 3 -GN have the same scaling dimensions at leading-order in 1/N : (3.9)

Real-space kernels
We now turn to the next-to-leading-order term in the free-energy expansion in Eq. (2.18). The free-energy correction is related to the second variation of the action by where r ≡ d 3 r g(r) and we defined the following kernels , (4.4) where S eff is defined in Eq. (2.12). The remaining scalar-gauge kernel 6 with mixed a µ (r) and σ (r ) partial derivatives is obtained by exchanging coordinates r, r in H q µ (r, r ), which has mixed σ(r) and a µ (r ) partial derivatives; thus we write H q µ (r , r) in Eq. (4.1). In terms of the fermions in the original system, the kernels are given by D q (r, r ) = ψ(r)ψ(r)ψ (r ) ψ (r ) s.p. , (4.5) where ψ is a single fermion flavor and the current is This can be re-expressed in terms of the fermionic Green's function G q (r, r ) = ψ(r)ψ (r ) s.p. and its hermitian conjugate G † q (r, r ) = − ψ (r ) ψ(r) s.p. . The Wick expansion of the kernels yields where the cyclicity of the trace was used. The remaining prefactor 2N in Eq. (4.1) is cancelled as the fluctuation fields are rescaled σ, a µ → σ/ √ 2N , a µ / √ 2N to control the expansion. The free-energy correction is then obtained by integrating the field fluctuations in Eq. (4.1). It is convenient to subtract the q = 0 correction F (1) 0 = 0, 7 therefore we write the general correction as where we define the matrix kernel . (4.13)

Fourier transform
To compute the determinant operator, the kernels are expanded in terms of harmonics. For the gauge-gauge kernels, the vector spherical harmonics are introduced: As suggested by the notation, the B mode has zero divergence: ∇ · a B m (n) = 0, and the E mode has zero curl: ∇ × a E m (n) = 0. It is also useful to introduce 4-dimensional eigenfunctions, , (4.17) 7 Since the scaling dimension of the identity operator vanishes, we requires less regularization procedures and automatically takes care of gauge fixing subtleties since the Fadeev-Popov ghost contribution is independent of the background flux 4πq in QED 3 [61].
where X ∈ {T, E, B}. In this basis, the matrix kernel M q (r, r ) can be expanded as where we directly work on S 2 × R (i.e., taking the limit β → ∞ now) and where All of the arguments of the functions appearing in the matrix are ω. Note that the scalar-gauge kernel is imaginary, hence the reason for the signs in the first column of the matrix. This last point is shown explicitly in App. B.
This kernel can be simplified by using CT invariance.
The auxiliary boson φ, a pseudo-scalar, and the E, T modes of the gauge field are antisymmetric under CT , while the B modes are symmetric under CT . This implies that the following kernels vanish 8 K q,T B (ω) = K q,EB (ω) = H q,B (ω) = 0. (4.20) The U(1) gauge invariance also enables the kernel to be simplified. Using the conservation of the U(1) current, ∇ µ J µ (r) = 0, in Eqs. (4.6-4.7), it follows that It should also be noted that among vector spherical harmonics, only a T µ, m (n) is defined for = 0. In this case, the only remaining gauge-gauge kernel is K q,T T 0 (ω), and it vanishes by gauge invariance. The computations to obtain these gauge invariance conditions are shown in App. C. Using all these simplifications, the monopole scaling dimension correction is given by where ω ≡ ∞ −∞ dω/ (2π) and we defined Note that H 0,T (ω) = 0 and thus it does not appear in the denominator. By turning off the GN interaction in Eq. (4.24), the scalar-scalar kernel D q (ω) and the scalar-gauge kernel H q,T (ω) do not contribute, and the monopole scaling dimension correction in QED 3 [16] is recovered: 8 This result was also checked explicitly with the same method giving the values of non-vanishing kernels.
One can alternatively deactivate the gauge field, keeping only the scalar-scalar kernels, and obtain the pure GN model. Despite the absence of a gauge field in this model, one can still introduce an external gauge field with the 4πq flux, define a correlation function on this background configuration, and obtain the related critical exponent. This was notably achieved for the O(N ) model in Ref. [62]. In a forthcoming publication, we shall also explore this avenue in the pure-GN model.
The relevant kernel Fourier coefficients to compute the monopole scaling dimensions in (4.24) and (4.27) are found by inverting Eq. (4.18): where the second coordinates are fixed to τ = 0 andn = z without loss of generality, and normalized coordinates a, a are introduced.

B. Anomalous dimensions
The anomalous dimensions of monopole operators (4. 24, 4.27) are computed in this section. To do so, the kernel coefficients in Eqs. (4.28, 4.30, 4.31) must be obtained. These coefficients are built with real-space kernels (4.9-4.11) that depend on the fermionic Green's function at the saddle-point. The Green's function is defined by the action of the Dirac operator on it: (4.32) The eigenkernels in Eqs. (4.28, 4.30, 4.31), involving the sums on spherical harmonics or vector spherical harmonics, will also be needed.

q = 0 kernels
We first compute the expressions in the denominator of the scaling-dimension corrections (4.24, 4.27), that is, the q = 0 kernel coefficients. The eigenkernel in the scalarscalar kernel (4.28) is just the sum of spherical harmonics, which is given by the addition theorem where cos γ ≡n ·n = cos θ cos θ + sin θ sin θ cos (φ − φ ) . (4.34) When working withn =ẑ, this is replaced by For the sums on vector spherical harmonics appearing in the gauge-gauge kernels (4.30, 4.31), a similar result is obtained in spherical coordinates a =θ,φ,τ in Eqs. (E.4, E.5) in App. E and is formulated with the same Legendre polynomial and its first and second derivatives. This reproduces a result from Ref. [63]. The real-space kernel for q = 0 is also needed. In this case, the Green's function takes a simple form which is simply the conformally transformed 3D flat space Green's function [16]: The real-space kernels can then be obtained in normalized spherical coordinates. Inserting this Green's function, along with the eigenkernels, in Eqs. (4.28, 4.30, 4.31), the resulting q = 0 kernel coefficients are (setting where integration by parts was used to eliminate derivatives of P (x). The differential operators acting on e iωτ P (x) can be replaced with the corresponding eigenvalues ∇ 2 S 2 → − ( + 1) and ∂ 2 τ → −ω 2 with further integration by parts. The remaining expressions contain Fourier transforms of the form r e iωτ P (x)X p which were obtained in the appendix of Ref. [63]. Using these results, the q = 0 kernel coefficients are simplified to (4.44) Note that we have reproduced the gauge-gauge coefficients K 0,E (ω) and K 0,B (ω) given in Ref. [16] by using the methods of Ref. [63].

Anomalous dimensions for q = 1/2
For the minimal magnetic charge, the eigenkernels in Eqs. (4.28-4.31) are formulated using the same expression (4.33, E.4, E.5) as in the last section. In particular, the gauge-gauge kernels will be worked out in normalized spherical coordinates. As for the real-space kernels (4.9-4.11), they depend on the q = 1/2 fermionic Green's function defined through Eq. (4.32). The spectral decomposition of the Green's function in terms of spinors with monopole harmonics components is shown in App. D 2. A generalized addition theorem for monopole harmonics involving the Jacobi Polynomials P (0,2q) (x) is then needed. Specifically, after taking the sum over the azimuthal quantum number, the Green's function for general q is given by [16] where the energies E q, were defined in Eq. (3.8) and where The phase e −i2qΘ comes from the generalized addition theorem and is defined in Eq. (D. 19), but it is not involved in the computation since it is always cancelled by the opposite phase of the Green's function hermitian conjugate. The Green's function can be inserted in Eqs. (4.9-4.11) to obtain the real-space kernels, which, along with the eigenkernels (4.33, E.4, E.5), are inserted in Eqs. (4.28-4.31) to compute the four kernel coefficients. Defining K q,D (ω) ≡ D q (ω) and K q,T (ω) ≡ H q,T (ω), the kernel coefficients K q,Z (ω) with Z ∈ {D, T, E, B} are given by where the prefactors are given by There is a sign error in the first term of the Green's function in Ref. [16] that we corrected here. This sign does not affect the conclusions in Ref. [16].
the integrals for scalar-scalar and scalar-gauge kernels are 52) and the integrals for gauge-gauge kernels reproduce the expressions obtained in Ref. [16] 10 These integrals can be performed exactly, see App. F for more details. In the end, these quantities depend only on the angular momenta: I Z 1 ( , , ) and I Z 2 ( , , ). For = = q, this computation requires more care since both energies vanish, and the integral over time leading to the prefactor in Eq. (4.47) instead yields a Dirac delta function δ(ω). However, for = = q = 1/2, the term in the bracket simply vanishes. When only one of and have their minimal value q = 1/2, there is a non-vanishing contribution to the anomalous dimension. In this case, only I Z 1 contributes, since the prefactor in front of I Z 2 in Eq. (4.47) vanishes. For q = 1/2, the contribution of zero modes in Eq. (4.47) vanishes with = ω = 0, otherwise it is given by (4.57) The remaining contribution consists in a sum on nonzero modes , ≥ 3/2. The summand depends on I Z 1 ( , , ) and I Z 2 ( , , ) which are formed of three-J symbols in , and (F.6-F.8). Thus, one of the sums, say on , can be viewed as finite. Then, after taking the sum on , the remaining summand tends to a constant for large (4.58) Thus, for kernels with a non-zero asymptotic constant, the sum on will be divergent. This is regularized with a zeta function regularization (4.59) The sum above is then finite and is computed numerically up to a cutoff c . The remainder is approximated with a large expansion of the summand −α Z + k . Each power in the expansion is summed analytically from = c + 1 to = ∞ with a zeta function. The coefficients c (ω) are found by doing the expansion for a few fixed values of and deducing the general dependence on . It turns out that only even powers of 1/ have non-vanishing coefficients c (ω). We obtained the expansion up to k = 18. With this remainder, we found that a cutoff c = 300 + 1/2 was sufficiently large to achieve the desired precision goals. The first few terms of the remainders for general q are shown in App. G.
After performing the sums in Eqs. (4.57, 4.59), the kernel coefficients in Eq. (4.47) are computed and inserted in Eq. (4.24) (or Eq. (4.27) for the case of QED 3 ). The kernel coefficients in the denominator of the logarithm of the monopole anomalous dimension were obtained analytically in Eqs. (4.41-4.43). The remaining sum on and integral on ω are computed up to a relativistic cutoff [16] ( + 1) + ω 2 ≤ L (L + 1) . (4.60) We obtained the anomalous dimension with a cutoff up to L max = 65. A function of 1/L is then fitted to extract the value of the anomalous dimension as the full sum and integration are taken with L → ∞. Fig. 1 shows a quartic function fitted with the data from L ∈ [L max − 10, L max ].
The anomalous dimension of a charge q = 1/2 monopole in QED 3 -GN extracted from this fit is ∆ This reproduces the result in Ref. [16] up to a difference of order 10 −4 .
While we extrapolated the result for L → ∞ with a quartic fit, based on a cutoff of L max = 65, varying the maximal relativistic cutoff can change the last digit in the result quoted above. For instance, ∆ 1/2,QED3 | Lmax=50 = −0.03815. In App. I, we show how we computed anomalous dimensions for various L max and used the trend as L max → ∞ to estimate the anomalous dimensions and their errors. The error we quote in what follows reflects the uncertainties related to the extrapolations and not the precision of our computation which yields relatively negligible errors.
Using this method, the anomalous dimension of q = 1/2 monopoles at next-to-leading order in the 1/N expansion in QED 3 -GN is given by (7). (4.61) The scaling dimension of q = 1/2 monopole operators in QED 3 -GN is then given by 2N × 0.26510 + 0.118911(7) + O N −1 . In QED 3 , the correction we found is With this estimated uncertainty of our result, it is clearer that there is a small discrepancy when comparing our result with the correction −0.0383 computed in Ref. [16]. Trying to replicate the method in Ref. [16], we used a cubic fit with data L ∈ [5,45] and obtained −0.03823, which is closer to −0.0383. with components given by the product of two monopole harmonics (D. 16-D.18). Consequently, the real-space kernels are formulated as products of four monopole harmonics. As for the eigenkernels appearing in D q (ω), H q,T (ω), they are already expressed as the product of two spherical harmonics (4.28, 4.29). Only the gaugegauge kernels then need a reformulation. In this case, a different basis U µ ,m (n), V µ ,m (n), W µ ,m (n) for the vector spherical harmonics can be introduced. These are eigenfunctions with respective total spin j = − 1, and + 1 [16]. Most importantly, the components of these harmonics are simply given by spherical harmonics (see App. E 2). We discuss the relation with the previous basis shortly.
Before doing so, we note that, in this new formulation, the kernel coefficients are expressed as the integral of a product of four monopole harmonics and two spherical harmonics. To be more precise, half of these functions are conjugate harmonics, but they can all be expressed as harmonics with the following relation Just as we did in Sec. IV B, the primed coordinates can be fixed as τ = 0 andn =ẑ without loss of generality. As a result, half of the six harmonics are eliminated This removes every sum on azimuthal quantum numbers, which greatly simplifies the computation. There remains an integral over three harmonics (4.65) The explicit expressions for the kernel coefficients involve the sum of many such integrals and are not reproduced here.
Returning to the change of basis, the U, V, W vector spherical harmonics in the j = sector can be related to the harmonics previously introduced in Eqs. (4.14-4.16) by The Fourier coefficients can also be transformed in this basis: (4.67) The matrix of eigenkernels keeps the same structure thanks to the block-diagonal form of the transformation matrix R. This is expected, as we could also argue that the kernels K q,U V (ω) = K q,W V (ω) = 0 vanish because of CT invariance, as we did for K q,T B (ω) = K q,EB (ω) = 0. The relevant relations are then The first relation is found by comparing the bottomright components in Eq. (4.67) and using the definition of K q,B (ω) (4.26), whereas the second relation is found by taking the trace of Eq. (4.67), using the first result in Eq. (4.68) and the definition of K q,E (ω) (4.25). The kernels K q,E (ω) and K q,B (ω) can then be replaced in the scaling-dimension corrections (4.24-4.27) by their formulation in the new basis. For general charge, the regularization of the kernels presented in Eqs. (4.58, 4.59) is still valid : Regulator terms −1/ (2π) and −1/ (4π) can be used respectively for the scalar-scalar and gauge-gauge kernels, while the scalar-gauge kernel does not require regularization. The contribution of the zero modes using this method is also very straightforward and algorithmic. However, there seems to be additional contributions coming from the combination of zero modes in both Green's function, = = q. As discussed previously, this contribution is proportional to a Dirac Delta function δ(ω) instead of the energy prefactor in Eq. (4.47). This contribution vanishes once integrated over ω (see App. H). Numerical sums on are obtained up to c = 200 + q, 11 and the remainder is computed analytically with an expansion up to 1/ 18 . In this case, the coefficients of the expansion also depend on the charge q and are found by fixing and q for a few values.
The anomalous dimension ∆ (1) q for each charge are computed with a relativistic cutoff L max = 35 + q , where q ≡ Round(q). Here, the convention is that halfintegers are rounded to even numbers, e.g. 1/2 = 0 and 3/2 = 2. The value of L max is modulated with the charge q to ensure that a regime with a tail-like behaviour, as observed in Fig. 1, is attained for larger charges. For the three minimal charges, we used a larger relativistic cutoff: We use L max = 65 for q = 1/2 (as in the last section) and L max = 46, 47 for q = 1, 3/2. 12 We found that the results are robust as L max is increased and more precise (see App. I). The results for L ∈ [L max − 6, L max ] are used to fit a quartic function in 1/L to extrapolate the anomalous dimensions ∆ As before, the uncertainty in the scaling dimension is estimated by varying L max and estimating the anomalous dimension as L max → ∞. More details are shown in 11 We use a smaller cutoff for the general charge q which is still sufficient for the precision needed and less computationally intensive. 12 In these cases, we use c = 300 + q. For q = 1/2, we fit with L ∈ [Lmax − 10, Lmax] (as in the last section) whereas we fit the q = 1, 3/2 cases with L ∈ [Lmax − 6, Lmax] (as other charges in this section). The number of points used for the fits is discussed further in App. I. Table I. Leading-order and next-to-leading order in 1/N contributions to monopole scaling dimensions in QED3, QED3-GN, and QED3-Z2GN models. The latter model is discussed in Sec. VII. The leading-order result is the same in all models. The scaling dimension in a given model is (7) 0.102846 (9)  The next-to-leading order term in 1/N decreases the scaling dimension of monopoles in QED 3 whereas it increases for QED 3 -GN. That is, quantum corrections help to stabilize the QED 3 -GN model and destabilize QED 3 . To understand the difference between both cases, it is useful to write the scaling dimension as is the contribution in the QED 3 -GN anomalous dimension (4.24) coming exclusively from the pseudo-scalar field. 13 Computing ∆ (1) q,GN in the same way as we did for ∆ (1) q,QED3-GN and ∆ (1) q,QED3 , we found this contribution is positive ∆ in the logarithm is positive. The numerator is explicitly positive. As for the denominator, we note that ∆ (1) q,GN is real, meaning that D q (ω) and the D 0 (ω) must have the same sign. The latter q = 0 kernel is positive, as seen from Eqs. (4.41, 4.44), meaning that D q (ω) > 0. The same goes for K q,E (ω), thus the denominator is positive D q (ω)K q,E (ω) > 0. This does not come as a surprise as these kernels are diagonal entries in the Hessian matrix developed around a minimum saddle-point. The scalar-gauge kernel thus gives a positive contribution to the anomalous dimension in QED 3 -GN. It then must be that the QED 3 -GN monopole anomalous dimension is positive, ∆ In the CP N −1 model, similar relations between the different contributions to the monopole anomalous dimension are found. A positive contribution coming only from the auxiliary boson was found numerically in Ref. [62]; the correction from the mixed scalar-gauge kernel can be deduced as positive [63; 64], and the total anomalous dimension of monopoles in CP N −1 numerically found in Ref. [64] is also positive.

Convexity conjecture
It was recently conjectured that CFT operators charged under a global U(1) symmetry respect the following convexity relation ∆((n 1 + n 2 )n 0 ) ≥ ∆(n 1 n 0 ) + ∆(n 2 n 0 ), (4.71) for some positive integer n 0 of order 1 [65]. We test this conjecture using the monopole operators that are charged under U(1) top . Here, n 0 , n 1 , n 2 are integers, where in our notation ∆(2q) ≡ ∆ q . Using the scaling dimensions we obtained in Tab. V and extrapolating to finite N , we find this relation is respected for the monopoles under consideration in QED 3 , QED 3 -GN (and also for the case QED 3 -Z 2 GN presented later on) for any 2N ∈ Z + starting from n 0 = 1, i.e. the minimal possible value.

V. LARGE-CHARGE UNIVERSALITY
In CFTs with a global U(1) symmetry, the related charge q can be used as an expansion parameter by using effective-field-theory methods. It was shown that the lowest scaling dimension among charge-q operators has the following expansion at q 1 [54]:  16 Since QED 3 and QED 3 -GN monopoles have the same leading-order scaling dimensions, as discussed in Sec. III, this also applies to QED 3 -GN monopoles.
Using the monopole anomalous dimensions ∆ q , the O q 0 coefficient γ can be computed. This was done for the CP N −1 model in Ref. [68], where ∆ (1) q was obtained for a hundred charges q = 1/2, 1, . . . , 50 and the expected expansion (5.1) is fitted numerically to extract γ. A similar computation is performed here for monopoles in the QED 3 -GN and QED 3 models. We fit all monopole anomalous dimensions in QED 3 -GN and QED 3 shown in Tab. V by using the fitting function in Eq. (5.1) with powers down to q −1 [66]. The fits and the anomalous dimensions are shown in Fig. 3; note that the errors in the values of the anomalous dimension are smaller than the dots in the figure. Including more powers in the fitting function would yield significantly larger errors in the estimation of γ.
The value of γ obtained for each theory is consistent 14 Here, the parameter N is used to designate either N complex boson flavors or 2N fermion flavors. 15 These models are also known as QED 3 −chiral XY GN and QED 3 −chiral Heisenberg GN models, respectively. 16 While Ref. [63] discusses only the O q 3/2 term of the large-q expansion, it is straightforward to use their analytical results to verify that no O N q 0 term is present.
with the expected universal value (5.2) This is a nice consistency check of the anomalous scaling dimensions obtained in the last section. The universal coefficient of the scaling dimension of U(1)-charged operators in Eq. (5.1) can also be formulated with the following sum rule [54] The RHS results from a cancellation of order O q 3/2 and O q 1/2 terms on the LHS. In comparison, the error on the LHS coming from the triplet ∆ q−1 , ∆ q , ∆ q+1 is comparatively large, even more so as the coefficients in front of the scaling dimensions, which are of order q 2 , become larger with increasing q. The resulting errors are too important to obtain a reliable fit of the large-q behaviour of this sum rule. Nevertheless, we did observe a good qualitative agreement for the scaling dimensions shown in Tab. V.

VI. CFT DUALITY: QED3-GN AND CP 1 MODELS
Another interesting application of our results concerns the duality between the QED 3 -GN model with 2N = 2 two-component Dirac fermion flavors and the CP N −1 model with N = 2 complex boson flavors [53]. Crucially, the duality between these models implies an emergent SO(5) symmetry. The following SO(5) multiplet in the is dual to the following multiplet in the CP 1 model  Table II. Operators in the SO(5) 5 multiplets (6.1, 6.2) and their scaling dimensions. "VBS" and "Néel" make reference to operators whose scaling dimensions are obtained numerically on lattices. The results for monopole operators are obtained by using the state-operator correspondence at next-to-leading order in 1/N . The scaling dimension of the auxiliary boson φ in QED3-GN was obtained at order 1/N using the mean of Padé and Padé-Borel [0/1] resummations (non-resummed scaling dimension are unphysical). The scaling dimension of the fermionic monopole operator can also be resummed to (0.59, 0.68), but not in the bosonic case. The operator z † σz designates any of the boson bilinears, i.e. flavor spin-1 in the bosonic side. It was obtained at order 1/N 2 in Ref. [69] and using functional renormalization group in Ref. [70]. within a multiplet should be equal, while operators identified by the duality should also have the same scaling dimension. Putting this together, this means that all the operators above should have the same scaling dimension. A decent agreement was already observed in Ref. [31], but the scaling dimension of QED 3 -GN monopoles were obtained only at leading-order in 1/N , with ∆ M f 1/2 = 0.53.
Updating the comparison with the next-to-leading order correction, we find ∆ M f 1/2 = 0.65, which gives an even better agreement. For instance, the scaling dimension of the q = 1/2 monopole on the CP 1 side obtained at nextto-leading in 1/N is given by ∆ M b 1/2 = 0.63 [63; 64]. In contrast, if we extrapolate the large-N QED 3 result to 2N = 2, we obtain a scaling dimension of 0.49, which is further from the CP 1 result, as expected since the two CFTs are not related by duality. The scaling dimensions of the other operators in the duality also show a good agreement coming from both analytical and numerical studies, as shown in Tab. II.
In the same way, monopoles with the second smallest charge q = 1 were argued to be part of the symmetric traceless 14 representation of SO(5) [44; 69]. The various relevant scaling dimensions obtained with analytical methods are also compared in Tab. III. Again there is a very good agreement between the scaling dimension of monopole operators, with ∆ M f 1 = 1.58 and ∆ M b 1 = 1.50. The agreement is weaker with other operators, but by taking into account Padé and Padé-Borel resummations the duality prediction seems quite reasonable. The scaling dimension related to auxiliary bosons ∆ φ 2 and ∆ λ obtained using the large-N have greater discrepancy with ∆ ∼ 1, but these expansions are not very well controlled. However, the same can be said about the monopole oper- Table III. Operators in the SO(5) symmetric traceless 14 multiplets and their scaling dimensions predicted to be equal according to the duality between QED3-GN |2N=2 and CP 1 models. The scaling dimensions presented are obtained analytically with the large-N expansion. Padé and Padé-Borel [0/1] resummations are shown in parenthesis (apart from Ref. [47], resummations are not obtained in the references cited). The symbol "×" indicates unphysical results, i.e., negative scaling dimensions. The operator λ is the Lagrange multiplier field on the CP 1 side. Results for monopole operators are obtained using state-operator correspondence at order N 0 , while other results were obtained at order N −1 . The resummed value for ∆ψ σψ was obtained in Ref. [47] and is the same for ∆ (z * σz)(z * σz) at this order.
Ref. ator on the bosonic side. Overall, the duality for the 14 representation of SO (5) is not as convincing as for the 5, but still reasonable for perturbative results. The scaling dimension of the Lagrange field obtained using the functional renormalization group also agrees reasonably well ∆ λ = 1.21 [70]. The situation becomes more puzzling when the critical exponents in Tab. III are compared to numerical lattice results. The apparent consistency observed in the analytical results (at least for operators that do not need resummation) does not hold for the numerical lattice results. Specifically, we compare the analytical results to the correlation length exponent ν obtained in many numerical studies of the CP 1 model. This exponent is related to the Lagrange field scaling dimension as ∆ λ = 3 − 1/ν. Its value has varied greatly among many numerical works. Earlier results indicate that ∆ λ ∈ [1.34, 1.67] [71; 72; 78; 79] which seems compatible with other scaling dimensions in Tab. III. However, unusual scaling behaviour and the "drifting" of ν with increasing lattice size [75] motivated further studies, and lower scaling dimensions have been found. The wide range of values obtained are shown in Tab. IV. Notably, a scaling dimension going down to ∆ λ = 0.80(1) by considering the presence of a second length scale [80].
The varying results among different lattice studies were also interpreted as a hint for a weakly first-order transition. This possibility has been discussed [69; 81] in a field theory context where the dual models QED 3 -GN | 2N =2 and CP N −1 | N =2 are possibly complex CFTs emerging from the collision of fixed points as the number of matter flavors is lowered below a critical level. On the other hand, our analytical analysis shows there is still consistency among scaling dimensions on both sides of the du- ality. This may imply that the duality can still give valuable information, even if the CFT is non-unitary.
A similar tension between the results from field theory and lattice models was observed in Ref. [83] where the QED 3 | 2N =2 model was studied using conformal bootstrap. The duality to the easy-plane CP 1 model conjectured in Ref. [53] implies a self-duality and an emergent O(4) symmetry on both sides. While the conformal bootstrap study of QED 3 | 2N =2 is consistent with the selfduality and the emergent symmetry, it contradicts results from the lattice study of the easy-plane CP 1 model [84].
An interesting approach to understand these discrepancies could be that of pseudo-criticality, that is a weakly first-order transition with a generically long correlation length. In Ref. [85], a Wess-Zumino-Witten model in 2+ dimensions, with target space S 3+ , with global symmetry SO(4 + ) has been shown to exhibit this behaviour and consistent with numerical results in the literature. A crucial point was that the physical dimension d = 3 is close to the critical dimension d = 2.77 where fixed points collide. Pseudo-criticality was also found in a loop model describing the easy-plane Néel-VBS transition [86].

Higher charge
The duality can be tested further by comparing monopoles on both sides of the duality. First, the relation between minimally charged monopoles is further discussed. This relation is simpler to see with the appropriate sub-models. A duality between QED 3 | 2N =2 and easy-plane CP 1 was formulated by including additional external gauge fields B µ , B µ and Chern-Simons terms [53; 87] where b µ and a µ are the dynamical gauge fields in bosonic and fermionic models, respectively. By inspecting the charges under the external gauge fields (q B , q B ), we can identify the following bosonic operators the scaling dimensions of higher charge monopoles can also be compared. The OPE of two q = 1/2 monopole operators yields the expansion over q = 1 operators where the ellipsis stands for other primary operators with larger scaling dimensions. By definition, M 1 (x) has the smallest scaling dimension in the q = 1 topological sector. We can then identify the scaling dimension of q = 1 monopole operators on both sides of the duality ∆ f q=1 = ∆ b q=1 . This is expected, as these monopoles are conjectured to be components dual SO(5) symmetric traceless 14 multiplet [44; 69]. Using the same logic for higher charge monopoles, we find more generally that Comparing our results for QED 3 -GN | 2N =2 monopoles in Tab. I to CP 1 monopoles in Ref. [68], we obtain a good agreement for higher charges, as shown in Fig. 4. For larger charges, the relative difference tends to 10%. This is a great improvement compared to the results obtained with only the leading-order scaling dimensions: the behaviour is similar and the asymptotic relative difference for large q is 76% instead.

VII. TRANSITION TO Z2 SPIN LIQUIDS
In this section, we consider quantum critical transitions to Z 2 spin liquids. The pairing of the spinons gaps out the gauge field through the Higgs mecanism. We begin with the most symmetric pairing interaction, and we then discuss the more general case where the pairing further breaks the flavor symmetry.

A. Symmetric Z2 spin liquid
The transition out of the U(1) DSL to a Z 2 spin liquid can also be studied with a gauged Gross-Neveu model [30; 49-51]. The Lagrangian describing this tran-sition is written in euclidean flat spacetime as where φ is a complex scalar that decouples a quartic superconducting pairing term for the fermions. The interaction term included preserves Lorentz invariance. As φ describes Cooper pairs, it transforms asψ i iγ 2ψ T i under U(1) gauge transformations. The Yukawa interaction term in the above equation is thus gauge invariant. In the Z 2 QSL, φ acquires an expectation value, which Higgses the gauge field, leading to a gapped s-wave superconducting state for the Dirac fermions.
We recall that the Dirac conjugate is defined by ψ i = ψ † i γ 0 and A q µ is the external gauge field that sources the flux of 4πq. The gauge-covariant derivative for the external gauge field A q µ on a curved spacetime is defined in Eq. (2.13). We now introduce the Nambu spinor defined as The transpose of the Nambu spinor is given by . Thus, the fermionic action can be expressed as where the (inverse) Nambu Green's function is As in Sec. IV, the fields φ and A are expanded about their saddle-point values as Here, G −1 0 is the bare inverse Green's function, determined from the gauge covariant derivative term involving A q µ in Eq. (7.4), and X σ and X a are given by Integrating out the fermions then gives the effective action as DΥ exp(−S Υ ) ≡ exp(−S eff ), where S eff = − 1 2 (2N )Tr log G −1 . We let Tr denote a "trace" over all relevant degrees of freedom, whereas tr denotes a trace over spinor components. To compute the effective action, we express the fermionic action as a quadratic form in the fluctuation fields and perform a Gaussian functional integral over a and σ. The linear terms in a and σ vanish due to the saddle-point conditions for A and φ. Thus, to quadratic order, the effective action becomes The fluctuations are O(1/ (2N )), thus they cancel the prefactor 2N . The second and third terms involve both the gauge field and the scalar field. By taking the trace over the Nambu matrix structure, these terms are found to vanish. Indeed, these terms must vanish from gauge invariance. Hence, only the gauge-gauge and scalar-scalar kernels contribute, which is in contrast to the QED 3 -GN case where mixing between the two sectors exists. After performing the trace over the Nambu indices, the scalarscalar kernel is where the scalar kernel is the same as in Eq. (4.9). Similarly, the gauge-gauge kernel is the same as in QED 3 : Tr G 0 (r , r)X a (r)G 0 (r, r )X a (r ) = 1 2 r,r a µ (r)K µν (r, r )a ν (r ), (7.9) where the gauge-gauge kernel is the same as in Eq. (4.10). Combining these two results we find that the fluctuation action (obtained after integrating out σ and a) is just the sum of twice the pure GN and the QED 3 results. The anomalous dimension for the minimal charge q = 1/2 is thus deduced to be (1) GN = 0.102846(9). (7.10) The value and the error are estimated in the same way as described in Sec. IV and App. I by finding an expression similar to Eq. (4.24) for the QED 3 -Z 2 GN case. The result is surprisingly close to the QED 3 -GN case, although the quantum fluctuations possess a different structure at the two transitions. In Table I, we give the answer for higher q. It can be seen that the values of the anomalous dimensions for q > 1/2 for the CSL and Z 2 QSL are not as close as in the case of the minimal charge. By using anomalous dimensions up to q = 13 shown in App. J, one can again confirm the value of the universal coefficient for CFTs with a U(1) global symmetry as described in Sec. V: In this case, we find γ QED3-Z2GN = 0.98(7) × γ U(1) .

B. More general Z2 spin liquids
Let us now consider a more general (superconducting) pairing interaction given by Here, C is the same as in the previous subsection -an antisymmetric and unitary matrix with indices denoting the Dirac indices of the spinor ψ. The additional term in the interaction M I represents a "flavor" matrix, which is symmetric and has indices in the valley and spin spaces. We shall consider some simple concrete examples of M later in this section. The index I, which is implicitly summed over, corresponds to the number of charged scalar fields, i.e., the competing pairing channels. In the previous section I = 1 and M corresponds to the identity operator.
The analysis for this more general pairing interaction follows the same lines as before. The Nambu spinor is the same as in Eq. (7.16), however, the inverse Green's function now becomes Here we suppose that M is hermitian. The inverse Green's function can be expanded again using the large-N formalism, and the effective action can be similarly expressed as in Eq. (7.7). The terms involving both the gauge field and the scalar field again vanish due to gauge invariance, and the gauge-gauge term is the same as in Eq. (7.9). The scalar-scalar kernel is now In the previous section, M I was the identity operator and so tr[(M I ) 2 ] = 2N , which leads to the previous result for the anomalous dimension in Eq. (7.10). As another example, consider the case where the pairing interaction is of the form φ x ψ T Cσ x ψ + φ z ψ T Cσ z ψ; note that we cannot use σ y since it is not symmetric. In this case numbers are expected to have different scaling dimensions, resulting in a hierarchy of monopoles [33]. As our present formalism selects only the monopole with the smallest scaling dimension, a constraint on the flavor quantum numbers would be needed to describe other monopoles. Moreover, since the pairing field cannot have an expectation value for gauge invariant monopoles, it is expected that next-to-leading corrections are necessary to observe this hierarchy. Generalizing Ref. [33] to quantify this effect would be an interesting avenue to explore. To understand these more general Z 2 spin liquids, we analyze the pairing Hamiltonian in further detail. In particular, here we shall focus on the Bogoliubov-de Gennes (BdG) Hamiltonian for the mean-field description of Eq. (7.11). In the preceding section we formulated the theory in terms of a Euclidean Lagrangian description. Since the Hamiltonian H is the time-component of an energy momentum tensor, it is necessarily a non-Lorentzinvariant entity. Thus, here we shall use Minkowski spacetime to perform the analysis, which enables standard field theory methods to determine H from L.
We define H = πψ − L, where π is the canonical momentum conjugate to ψ. From Eqs. (7.1) (without the gauge field) and (7.11), we construct the following Hamiltonian Our choice of gamma matrices is given by γ µ = (τ z , iτ x , iτ y ); this definition is consistent with the Clifford algebra with a mostly minus metric. Here, C = iτ y . We define the Nambu spinor χ(k) by .

(7.16)
In terms of the Nambu spinor χ(k), the Hamiltonian has the following form in momentum space: In general, the Hamiltonian can be expressed as where H(k) is the BdG matrix. In the simplest case we have 2N = 2, that is, there are 2 spin degrees of freedom, and the BdG Hamiltonian is 8 × 8. In the previous section we considered the case where the pairing matrix is the identity, M I = ∆σ 0 , where ∆ is the finite value of the pairing term in the Z 2 spin liquid phase.
Here we contemplate some simple examples of pairing matrices, namely M I = ∆σ x or M I = ∆σ z . For these classes of spin liquids, the BdG matrix is of the form Here σ l = σ 0 , σ x , or σ z . The various matrices appearing above correspond to the spin indices, Nambu indices, and finally the Dirac indices, respectively. For the simple case where σ l = σ 0 , the eigenvalues are given by E(k) = ± k 2 + 4 |∆| 2 , with a fourfold degeneracy.
These eigenvalues are exactly the same as in the Fu-Kane model at half filling [88]. Indeed, in the Fu-Kane model the Hamiltonian is a 4 × 4 matrix with spin and Nambu indices. Here we have two copies of this model Hamiltonian for each species of spin. In the case where σ l = σ x or σ z , we also have the same dispersion. In general, this dispersion describes a gapped Z 2 spin liquid. The QCP in the present model is generally thought to be well defined at modest values of N ; we can incorporate N copies of the valley degrees of freedom and obtain the same (copies) of the eigenvalues.
For the case where we have two competing channels, with twofold degeneracy. In the case where either ∆ x or ∆ z is equal to zero, we recover the previous result for the eigenvalues. Interestingly, for the specific case where |∆ x | = |∆ z | and arg(∆ x ) − arg(∆ z ) = ±π/2 (mod π), there are gapless Dirac cones, ±k. For such pairings, we thus have a gapless Z 2 spin liquid with massless relativistic fermions, and a gapped Z 2 gauge field.
This last case can be reformulated as a pairing term given by M I = |∆| d · σ y σ with d = e iϕ (1, 0, ±i). This expression is similar to the time-reversal-breaking triplet state defined with d = (1, i, 0), notably used to describe an LaNiC 2 compound [89], as well as a potential order parameter in twisted bilayer graphene [90]. This situation occurs when the general condensates ∆ 0 , ∆ z in M I = ∆ 0 σ 0 + ∆ z σ z are aligned in the complex plane, arg(∆ 0 ) − arg(∆ z ) ∈ {0, π}.
An important quantity for gapped systems is the Chern number [91], which corresponds to the flux of the Berry curvature in the Brillouin zone (BZ). In the presence of a gap the Chern number is an integer and describes a topological property of the system. Non-zero Chern number indicates broken time reversal symmetry, however, the converse is not always true. For a 2D system, it is defined by Here, B is the Berry curvature defined by B(k) = i (∇ k × A(k)) z , where A = n u n (k)|∇ k u n (k) with u n (k) a normalized eigenstate of the Hamiltonian, and the sum running over the occupied bands. Here we are considering a continuum theory where the bandstructure has a power-law dependence on momentum. In a lattice calculation, where the bandstructure is defined in the first BZ and involves trigonometric functions, the Chern number is well defined. To incorporate such physics, we extend the continuum model to a simple lattice dispersion where we replace k x , k y with sin(k x ), sin(k y ). This procedure does necessarily add additional Dirac points in the BZ. For the points in parameter space where there is a gap, we find that the Chern number is well defined and at all such points we obtain C = 0. This is expected when the system has time-reversal symmetry, which happens when both condensates ∆ x and ∆ z are imaginary. Other points that break TRS are connected without gap closing and thus are also expected to have a vanishing Chern number.

VIII. OTHER PHASE TRANSITIONS
The QED 3 -U(N ) × U(N ) GN, QED 3 -chiral XY GN, and QED 3 -cHGN transitions are also described with GN models, where the fermionic quartic interaction is respectively decoupled with N b = 1, 2, 3 real auxiliary bosons where the sum over 1 ≤ I ≤ N b is implicit, and µ I are Pauli matrices acting on a two-dimensional flavor subspace where As mentioned previously, the QED 3 -cHGN model describes the transition from a U(1) DSL to an AFM on the kagome lattice. In this case, the Pauli matrices µ act on magnetic spin subspace. As for the chiral XY interaction, taking µ x , µ y to act on a valley subspace, this describes the transition to a VBS order parameter. These transitions were observed for Monte-Carlo simulations on a square lattice where by tuning gauge-field fluctuations, the U(1) DSL is driven to either an AFM or VBS order, depending on the number of fermion flavors [92; 93]. A theoretical study that elucidated the field theory for the transition to an AFM was performed in Ref. [94] (see also Refs. [28; 31] for earlier studies of this model), while the field theory for the transition to the VBS was outlined in Refs. [95][96][97]. In this work, we used the appellation QED 3 -GN to designate the model with U(2N ) symmetry, following the convention of Refs. [47; 53], notably. The Pauli matrix in this case acts on valley subspace. However, the label was also used in the literature to refer to the U(N ) × U(N ) symmetric model, see Ref. [45] for instance. Both variations of QED 3 were considered in Refs. [48; 69]. As shown in Ref. [31], the auxiliary boson in these cases has a non-vanishing expectation value: |φ| = 0 in the monopole background on S 2 × R. This is also true for other choices of Pauli matrices, not only the specific one prescribed before to describe specific universality classes. For instance, the µ considered for the QED 3 -cHGN universality class could also act on valley subspace, in which case the order parameter is odd under time reversal.
There is still a non-vanishing expectation value of the auxiliary boson. Consequently, the Green's function in Eq. (4.45) must be modified to include a non-zero mass for the fermions. Using the addition theorems for spinor monopole harmonics needed to compute the zero-mass Green's function should be sufficient for this adaptation. Real-space kernels like in Eqs. (4.9-4.11) would also include Pauli matrices for traces on the magnetic spin subspace, with the number of kernels to compute increasing accordingly with the number of auxiliary bosons N b .

IX. CONCLUSION
We obtained the scaling dimension of monopole operators at the QCP between a U(1) DSL and two types of topological spin liquids, namely the CSL and a general class of Z 2 QSLs, at next-to-leading order in a 1/N expansion. The most relevant monopole operator in the CSL case has a minimal charge q = 1/2 and a scaling dimension ∆ 1/2,QED3-GN = 2N × 0.26510 + 0.118911 (7), while the analog scaling dimension in the case of the simplest Z 2 QSL is ∆ 1/2,QED3-Z2GN = 2N × 0.26510 + 0.102846 (9). For the other general Z 2 spin liquids, where the spin and valley flavour interaction was included, we obtained a general expression for the anomalous dimension. Since the spin/valley interaction reduces the size of the flavour group, an interesting question is what type of hierarchy the monopoles will have and how can one observe this in the monopole scaling dimensions. We also rederived the QED 3 monopole scaling dimensions and found small discrepancies, e.g., the q = 1/2 anomalous dimension is −0.038138(5) instead of −0.0383 and so on for other charges up to q = 5/2 [16; 61]. This will also lead to corrections to the anomalous dimensions of certain monopole operators in QCD 3 with non-abelian gauge groups, such as U(N c ), where the QED 3 anomalous dimension makes its appearance [61].
With these anomalous dimensions, we obtained a fit in the topological charge q and compared the O(q 0 ) coefficient with the universal value obtained in a large-charge expansion for operators charged under a global U(1) symmetry [54]. We obtain the expected value γ = −0.0937 in QED 3 , QED 3 -GN and QED 3 -Z 2 GN. We also revisited the conjectured duality between QED 3 -GN | 2N =2 and in CP 1 models [53]. Notably, the q = 1/2 monopole scaling dimensions in QED 3 -GN | 2N =2 agree very well with the scaling dimensions of other operators that are predicted to be equal under the duality. Specifically, the anomalous dimension obtained in this work greatly improves this agreement. We also argued that all monopoles with equal charges should have the same scaling dimensions in the QED 3 -GN | 2N =2 and CP 1 models. Using next-to-leading order results for both models, we obtain an agreement that is better for a minimally charged monopole, with a relative difference of 3%. As the topological charge increases, this difference increases and eventually saturates at 10% for q → ∞.
It would be interesting to study monopole operators in the other gauged GN models that we briefly discussed, notably the model describing the transition to an AFM [31][32][33]. Another interesting aspect to consider that was not included in this work is the case of monopole operators in the pure-GN model. It is a straightforward adaptation to write out the monopole anomalous dimensions in this case and use the results of this work to obtain them. Although there is no U(1) top due to the absence of a gauge field in this model, these objects still have useful applications. This notably motivated the study of monopoles in the bosonic O(N ) model [62; 98]. A study of the GN global monopoles and some of their applica-tions will appear in a forthcoming work.
Taking the divergence of the kernels, we obtain where ω ≡ dω/(2π). Requiring gauge invariance and setting these divergences to 0, we obtain the following relations Let us check one important relation following from gauge invariance: This is easily verified for q = 0 where closed forms of the kernels are easily obtained as discussed in Sec. IV B 1. This is a bit more involved when q = 0. We have expressions where the dependence on ω is easily isolated, taking the following form In the RHS of the gauge invariance condition in Eq. (C.12), we may reexpress the ω-dependent function as (C.14) Upon integration over ω, the first term is a simple divergence that can be regularized away. In the presence of test function f (ω), this contribution is dωf (ω) × 1 and it again vanishes provided the test function is convergent with no poles. The gauge-invariance condition in Eq. (C.12) can then be written by simply comparing the finite parts This last relation following from gauge invariance was verified for q = 1/2. Note that gauge invariance also implies that = 0, (C. 16) which is also verified by direct computation.
Appendix D: Green's function

Eigenvalues of determinant operator
In a general basis, the gauge covariant derivative acting on a spin-1/2 spinor on spacetime M will take the form where Ω µ is the spin connection transporting the fermion fields on spacetime M. On a flat spacetime M = R 3 , there are also spin connections in spherical coordinates that can be eliminated with a unitary transformation [99]. In this case, the covariant derivate ∇ µ=r,θ,φ can be traded for a normal derivative ∂ µ=r,θ,φ / D Proceeding with the Weyl transformation ψ → e −τ ψ, g µν → e −2τ discussed in Eq. (2.6), the Dirac operator on S 2 × R is given by [15] / D To diagonalize this operator, we introduce spinor monopole harmonics These spinors diagonalize the following generalized total spin and angular momentum operators J 2 q , J z q , L 2 q . In particular, the spinor monopole harmonics S ± q, ,m have a total spin j = ± 1/2. In the j = − 1/2 basis, the Dirac operator mixes the two types of spinors and simply becomes a matrix with c-number entries [15]: Here, the τ i are the Pauli matrices acting in the j = − 1/2 basis, i.e., they mix the components S + q, −1,m and S − q, ,m . For = q (we suppose a positive magnetic charge q > 0), only S − q,q,m exists and corresponds to a zero mode of the Dirac operator. By diagonalizing this matrix, we retrieve the eigenvalues used in Sec. III.

Green's function
The Green's function can be obtained with the spectral decomposition where ψ λ (r) are eigenspinors of the Dirac operator i / D A q ψ λ = E λ ψ λ forming a complete basis λ ψ λ (r)ψ † λ (r ) = δ(r − r ). With this formulation, the Green's function respects its defining equation of motion (4.32). We can simply keep working in the spinor monopole harmonics basis instead of further diagonalizing. The spectral decomposition of the Green's function in this basis is then The eigenvalue matrix is block diagonal, separating each j = − 1/2 sectors. To obtain a Green's function which has two particle-hole indices but which is a scalar with respect to the +/− structure described above, we take the left eigenspinor as a row vector in the +/− space, e −iωτ S + q, −1,m (n), e −iωτ S − q, ,m (n) . The action of the Dirac operator on the left eigenspinor is then given by where we used that N T q, = N q, , M T q, = M q, and M q, N q, = −N q, M q, . We can read the eigenvalue matrix from this relation and write Green's function 17 We can note that N −1 q, = N q, as this matrix squares to identity N 2 q, = −2 q 2 + E 2 q; τ 0 = τ 0 . Also, by noting that (ω + iM q, ). Then, the inverse matrix in the spectral decomposition becomes (ω + iM q, ) N q, (D.12) where we used that M q, N q, = −iE q; τ y . The spectral decomposition of the Green's function then becomes 14) The contour integral on ω is obtained with the residue theorem .
(D. 15) The spectral decomposition after the ω integration becomes By inserting Eq. (D.4), we obtain a 2 × 2 matrix whose components are pairs of monopole harmonics .
(D.18) 17 The inverse eigenvalue matrix here is N q, ω − iM q, −1 instead of N q, ω + iM q, −1 as described in Ref. [16]. This explains the benign different sign in our Green's function in what follows.
This formulation is used in Sec. IV B 3. For the minimal charge case, the Green's function can be further simplified by taking the sum on the azimutal quantum number [16], which yields Eq. (4.45) in the main text 18 . The phase appearing in Eq. (4.45) is given by [16] e −iΘ cos γ 2 = cos θ 2 cos Using the addition formula in Eq. (4.33), we obtain the results in Ref. [63].
This may be further simplified by using the differential equation for Legendre polynomials (E.6)

Second basis
By working in helical coordinates, the vector spherical harmonics U a m (n), V a m (n), W a m (n) introduced in Sec. IV B 3 take the following form Using a transformation matrix these harmonics can be rotated to cartesian coordinates e µ a Z a ,m (n), Z = U, V, W, (E. 12) which corresponds to the harmonics used in the main text.

Generalization
The kernel coefficients, with the "zero-zero mode contribution" explicitly included, can be written as In fact, it turns out that C D = −iC HT = −C KE , but this is not necessary for the argument that follows. Again, let us consider the calculation of the scaling dimension in QED 3 -GN near ω = 0. Once again, we can ignore the denominator: (2 + 1) ln K q,B (ω) D q (ω)K q,E (ω) + 1 + ω 2 ( + 1) F q,T (ω) 2 . (H.14) The K q,B (ω) is also regular and can be removed. Also, we already showed that − dω 2π ln D q (ω) = 0, we can use this to eliminate the = 0 contribution. We are left with (2 + 1) ln D q (ω)K q,E (ω) + 1 + ω 2 ( + 1) F q,T (ω) 2 .
(H. 15) Let us now consider the argument of the logarithm D q (ω)K q,E (ω) + 1 + ω 2 ( + 1) F q,T (ω) 2 = [C D reg KE + C D reg KE + 2 Re (C F T reg * F T )] δ(ω) (H. 16) Then, the scaling dimension correction near ω = 0 can be written as 4. Extrapolate the behaviour as L max → ∞ with a linear fit in 1/L max ,c 0 +c 1 /L max , using the five anomalous dimensions obtained. The fit and the anomalous dimensions are shown in Fig. 5 for the case q = 1/2. Additional points down to L max = 45 are also displayed in Fig. 5 to discuss the behaviour later on.
5. Compare the anomalous dimension obtained with the maximal value of L max , that we note in what follows as L * max (for q = 1/2, L * max = 65) and the extrapolated value at L max → ∞ and estimate the anomalous dimension as This result with the error bar is shown in Fig. 5 for the case q = 1/2 and is given by: 1/2,QED3 = −0.038138(5), ∆ 1/2,QED3-GN = 0.118911 (7).

(I.2)
We emphasize that the data used for the extrapolation in step 4 was itself the result of the extrapolation in step 2 (and step 3). The extrapolated value at L max → ∞ is therefore used as a guiding value. Had we taken a dataset with smaller values of L max , the fitted line would simply overshoot the one currently presented in Fig. 5. The anomalous dimension in Eq. I.1 would have more extreme values and thus a greater error.
To further characterize the effect that the size of L max has on the anomalous dimensions, we also consider the cases q = 1, 3/2 and q = 25/2 where we used a cutoff L max = 45 + q . 22 For those charges obtained with a larger cutoff, we can restrain our dataset and obtain the anomalous dimensions with L max = 35 + q . The results with L max = 45 + q are slightly more precise and very similar to those with L max = 35 + q . As shown in Fig. 6, the   drift of the anomalous dimension as L max is increased is very small relatively to the estimated errors, which indicates the stability of our method.
The same procedure was also used for different fitting functions k i=0 c i;k;Lmax L −i with higher polynomial order k = 5, 6, as shown in Fig. 7. We find a similar behaviour and more precise results. However, these fits demand a larger dataset. For larger q (and therefore larger maximal cutoff since the cutoff increases with q ), a similar behaviour remains. However, the size of datasets needs to be increased. This was also observed for q = 1/2 when comparing relativistic cutoffs L max of different sizes. It may indicate that errors are overfitted for smaller datasets with higher-order fits, as the effect is less important for k = 4. We used the quartic fit for all of the anomalous dimensions quoted in this work.    Appendix J: Monopole scaling dimensions for 1/2 ≤ q ≤ 13 Table V. Scaling dimension of monopole operators at leading-order and next-to-leading order in 1/N in QED3, QED3-GN and QED3-Z2GN models. The leading-order result is the same in all models. The scaling dimension is 2N ∆