Critical properties of the N\'eel/algebraic-spin-liquid transition

TL;DR

This paper investigates the quantum phase transition from Ne9el order to an algebraic spin liquid in 2D systems, identifying its universality class and calculating critical exponents using advanced field-theoretic methods.

Contribution

It establishes the universality class of the Ne9el/algebraic spin-liquid transition as chiral Heisenberg QED3-Gross-Neveu-Yukawa and computes critical exponents via epsilon and large-Nf expansions.

Findings

The transition is continuous with specific critical exponents.

The universality class is identified as chiral Heisenberg QED3-Gross-Neveu-Yukawa.

Critical exponents for thermodynamic and susceptibility quantities are calculated.

Abstract

The algebraic spin liquid is a long-sought-after phase of matter characterized by the absence of quasiparticle excitations, a low-energy description in terms of emergent Dirac fermions and gauge fields interacting according to (2+1)D quantum electrodynamics (QED$_3$), and power-law correlations with universal exponents. The prototypical algebraic spin liquid is the Affleck-Marston $\pi$-flux phase, originally proposed as a possible ground state of the spin-1/2 Heisenberg model on the 2D square lattice. While the latter model is now known to order antiferromagnetically at zero temperature, recent sign-problem-free quantum Monte Carlo simulations of spin-1/2 fermions coupled to a compact U(1) gauge field on the square lattice have shown that quantum fluctuations can destroy N\'eel order and drive a direct quantum phase transition to the $\pi$-flux phase. We show this transition is in the universality class of the chiral Heisenberg QED$_3$-Gross-Neveu-Yukawa model with a single SU(2) doublet of four-component Dirac fermions (i.e., $N_f=1$), pointing out important differences with the corresponding putative transition on the kagome lattice. Using an $\epsilon$ expansion below four spacetime dimensions to four-loop order, and a large-$N_f$ expansion up to second order, we show the transition is continuous and compute various thermodynamic and susceptibility critical exponents at this transition, setting the stage for future numerical determinations of these quantities. As a byproduct of our analysis, we also obtain charge-density-wave and valence-bond-solid susceptibility exponents at the semimetal-N\'eel transition for interacting fermions on the honeycomb lattice.

Figures

Tables

Table 1. Table 1: Padé and Padé-Borel resummations of the ϵ italic-ϵ \epsilon -expansion expressions for d = 3 𝑑 3 d=3 and N f = 1 subscript 𝑁 𝑓 1 N_{f}=1 . Approximants which are either singular in the domain 0 < ϵ < 1 0 italic-ϵ 1 0<\epsilon<1 , undefined, or negative, are denoted by × \times .

$N_f=1$	${\eta }_{\varphi }$	$1/\nu$	$\omega$	${\Delta }_{\overline{\Psi }\Psi }$	${\Delta }_{i\overline{\Psi }{\Gamma }_5\Psi }$
$\left[0/2\right]$	$\times$	0.334	$\times$	1.67	0.723
${\left[0/2\right]}_{\text{PB}}$	$\times$	0.636	$\times$	2.34	1.34
$\left[1/1\right]$	11.1	$\times$	0.526	3.17	$\times$
${\left[1/1\right]}_{\text{PB}}$	$\times$	$\times$	0.578	3.23	$\times$
$\left[0/3\right]$	$\times$	0.0826	$\times$	1.11	0.784
${\left[0/3\right]}_{\text{PB}}$	$\times$	0.534	$\times$	$\times$	1.20
$\left[1/2\right]$	39.4	$\times$	$\times$	0.393	0.776
${\left[1/2\right]}_{\text{PB}}$	$\times$	$\times$	$\times$	2.39	0.787
$\left[2/1\right]$	135	$\times$	0.894	$\times$	$\times$
${\left[2/1\right]}_{\text{PB}}$	$\times$	$\times$	0.862	$\times$	$\times$
$\left[0/4\right]$	$\times$	$\times$	$\times$	0.452	0.0986
${\left[0/4\right]}_{\text{PB}}$	$\times$	0.492	$\times$	2.05	1.10
$\left[1/3\right]$	1.09	0.208	0.0596	$\times$	0.724
${\left[1/3\right]}_{\text{PB}}$	2.26	0.236	0.664	2.34	$\times$
$\left[2/2\right]$	11.7	$\times$	1.71	$\times$	$\times$
${\left[2/2\right]}_{\text{PB}}$	$\times$	$\times$	$\times$	$\times$	$\times$
$\left[3/1\right]$	5.22	$\times$	1.27	$\times$	$\times$
${\left[3/1\right]}_{\text{PB}}$	5.29	$\times$	1.48	$\times$	$\times$

Table 2. Table 2: Padé and Padé-Borel resummations of the large- N f subscript 𝑁 𝑓 N_{f} expressions for d = 3 𝑑 3 d=3 and N f = 1 subscript 𝑁 𝑓 1 N_{f}=1 . Approximants which are either singular in the domain N f ≥ 1 subscript 𝑁 𝑓 1 N_{f}\geq 1 , undefined, or negative, are denoted by × \times . The exponents that are unknown beyond 𝒪 ​ ( 1 / N f ) 𝒪 1 subscript 𝑁 𝑓 \mathcal{O}(1/N_{f}) are denoted by − - , for these quantities only one approximant can be used.

$N_f=1$	${\eta }_{\varphi }$	$1/\nu$	${\Delta }_{\overline{\Psi }\Psi }$	${\Delta }_{i\overline{\Psi }{\Gamma }_5\Psi }$	${\Delta }_{i\overline{\Psi }{\Gamma }_3{\Gamma }_5\Psi }$
$\left[0/1\right]$	$\times$	0.4022	1.873	1.458	7.783
${\left[0/1\right]}_{\text{PB}}$	$\times$	0.5192	1.880	1.540	$\times$
$\left[0/2\right]$	$\times$	–	1.891	1.138	–
${\left[0/2\right]}_{\text{PB}}$	$\times$	–	$\times$	1.351	–
$\left[1/1\right]$	$\times$	–	1.889	$\times$	–
${\left[1/1\right]}_{\text{PB}}$	$\times$	–	1.887	$\times$	–

Table 3. Table 3: Padé and Padé-Borel resummations of the ϵ italic-ϵ \epsilon -expansion expressions for d = 3 𝑑 3 d=3 and N f = 1 subscript 𝑁 𝑓 1 N_{f}=1 in the chiral Heisenberg GNY model. Approximants which are either singular in the domain 0 < ϵ < 1 0 italic-ϵ 1 0<\epsilon<1 , undefined, or negative, are denoted by × \times .

$N_f=1$	${\Delta }_{\overline{\Psi }\Psi }$	${\Delta }_{i\overline{\Psi }{\Gamma }_5\Psi }$
$\left[0/2\right]$	2.55	1.79
${\left[0/2\right]}_{\text{PB}}$	2.66	2.01
$\left[1/1\right]$	$\times$	1.51
${\left[1/1\right]}_{\text{PB}}$	$\times$	$\times$
$\left[0/3\right]$	2.69	1.66
${\left[0/3\right]}_{\text{PB}}$	$\times$	1.92
$\left[1/2\right]$	2.65	1.55
${\left[1/2\right]}_{\text{PB}}$	2.69	1.56
$\left[2/1\right]$	2.62	1.54
${\left[2/1\right]}_{\text{PB}}$	2.61	1.54
$\left[0/4\right]$	2.02	1.57
${\left[0/4\right]}_{\text{PB}}$	2.54	1.87
$\left[1/3\right]$	2.57	1.24
${\left[1/3\right]}_{\text{PB}}$	2.67	$\times$
$\left[2/2\right]$	$\times$	1.52
${\left[2/2\right]}_{\text{PB}}$	$\times$	$\times$
$\left[3/1\right]$	2.51	1.53
${\left[3/1\right]}_{\text{PB}}$	2.52	1.53

Table 4. Table 4: Padé and Padé-Borel resummations of the large- N f subscript 𝑁 𝑓 N_{f} expressions for d = 3 𝑑 3 d=3 and N f = 1 subscript 𝑁 𝑓 1 N_{f}=1 in the chiral Heisenberg GNY model. Approximants which are either singular in the domain N f ≥ 1 subscript 𝑁 𝑓 1 N_{f}\geq 1 , undefined, or negative, are denoted by × \times .

$N_f=1$	${\Delta }_{\overline{\Psi }\Psi }$	${\Delta }_{i\overline{\Psi }{\Gamma }_5\Psi }$
$\left[0/1\right]$	2.508	1.816
${\left[0/1\right]}_{\text{PB}}$	$\times$	1.829
$\left[0/2\right]$	2.239	1.715
${\left[0/2\right]}_{\text{PB}}$	2.168	1.752
$\left[1/1\right]$	2.319	1.559
${\left[1/1\right]}_{\text{PB}}$	2.325	$\times$

Table 5. Table 5: PSG transformation properties of spin-triplet Dirac fermion bilinears.

	$\overline{\psi }𝝈\psi$	$\overline{\psi }{\tau }_1𝝈\psi$	$\overline{\psi }{\tau }_2𝝈\psi$	$\overline{\psi }{\tau }_3𝝈\psi$
$C_4$	$+$	$-\overline{\psi }{\tau }_2𝝈\psi$	$\overline{\psi }{\tau }_1𝝈\psi$	$+$
$R_x$	$-$	$+$	$-$	$+$
$T_x$	$+$	$+$	$-$	$-$
$T_y$	$+$	$-$	$+$	$-$
$𝒯$	$+$	$-$	$-$	$-$

Table 6. Table 6: PSG transformation properties of spin-singlet Dirac fermion bilinears.

	$\overline{\psi }\psi$	$\overline{\psi }{\tau }_1\psi$	$\overline{\psi }{\tau }_2\psi$	$\overline{\psi }{\tau }_3\psi$
$C_4$	$+$	$-\overline{\psi }{\tau }_2\psi$	$\overline{\psi }{\tau }_1\psi$	$+$
$R_x$	$-$	$+$	$-$	$+$
$T_x$	$+$	$+$	$-$	$-$
$T_y$	$+$	$-$	$+$	$-$
$𝒯$	$-$	$+$	$+$	$+$

Equations

$$H$$

$$+K\sum _{\square }\cos (\Delta \times \mathbf{\theta }).$$

$$H_0=\sum _{\sigma }{\int }_{\text{1BZ}}\frac{d^{2}k}{(2\pi {)}^2}{c}_{\mathbf{k}\sigma }^{†}h(\mathbf{k})c_{\mathbf{k}\sigma }^{},$$

$$\displaystyle h({\bm{k}})=-t\left(\begin{array}[]{cc}0&f({\bm{k}})\\ f^{*}({\bm{k}})&0\end{array}\right),$$

$$f(\mathbf{k})=1-e^{i(k_x-k_y)}+e^{-i(k_x+k_y)}+e^{-2i{k}_y}.$$

$$h(\pm \mathbf{Q}+\mathbf{p})=v_F(\pm {\mu }_1{p}_x+{\mu }_2{p}_y)+\mathcal{O}({\mathbf{p}}^2),$$

$$H_0\approx v_F\sum _{\alpha ,\sigma }\int \frac{d^{2}p}{(2\pi {)}^2}{\psi }_{\alpha \sigma }^{†}(\mathbf{p})\left({\mu }_1{p}_x+{\mu }_2{p}_y\right){\psi }_{\alpha \sigma }^{}(\mathbf{p}),$$

$$\displaystyle\psi_{+\sigma}({\bm{p}})=\left(\begin{array}[]{c}c_{{\bm{Q}}+{\bm{p}},A,\sigma}\\ c_{{\bm{Q}}+{\bm{p}},B,\sigma}\end{array}\right),\hskip 8.53581pt\psi_{-\sigma}({\bm{p}})=\left(\begin{array}[]{c}c_{-{\bm{Q}}+{\bm{p}},B,\sigma}\\ -c_{-{\bm{Q}}+{\bm{p}},A,\sigma}\end{array}\right).$$

$${\mathcal{L}}_0=\sum _{\alpha ,\sigma }{\overline{\psi }}_{\alpha \sigma }{\gamma }_{\mu }{\partial }_{\mu }{\psi }_{\alpha \sigma },$$

$$\mathbf{N}\equiv \frac{1}{\mathcal{N}}\sum _{\mathbf{R}}\left({\mathbf{S}}_{\mathbf{R}A}-{\mathbf{S}}_{\mathbf{R}B}\right),$$

$$\displaystyle{\bm{N}}\sim\sum_{\sigma\sigma^{\prime}}\int\frac{d^{2}p}{(2\pi)^{2}}\Psi_{\sigma}^{\dagger}({\bm{p}})\left(\begin{array}[]{cc}\gamma_{0}&0\\ 0&-\gamma_{0}\end{array}\right)\bm{\sigma}_{\sigma\sigma^{\prime}}\Psi_{\sigma^{\prime}}({\bm{p}}),$$

$$\displaystyle\Psi_{\sigma}=\left(\begin{array}[]{c}\psi_{+\sigma}\\ \psi_{-\sigma}\end{array}\right),\hskip 14.22636pt\sigma=\uparrow,\downarrow.$$

$$\displaystyle\Gamma_{\mu}=\left(\begin{array}[]{cc}\gamma_{\mu}&0\\ 0&-\gamma_{\mu}\end{array}\right),\hskip 14.22636pt\mu=0,1,2,$$

$$\mathbf{N}\sim \int d^{2}x\overline{\Psi }\mathbf{\sigma }\Psi ,$$

$${\mathcal{L}}_{2+1}=\sum _{\alpha ,\sigma }{\overline{\psi }}_{\alpha \sigma }{\gamma }_{\mu }{D}_{\mu }{\psi }_{\alpha \sigma }+\frac{1}{4}{F}_{\mu \nu }^2+g\mathbf{\varphi }\cdot \overline{\Psi }\frac{\mathbf{\sigma }}{2}\Psi +{\mathcal{L}}_{\mathbf{\varphi }},$$

$$\displaystyle\Gamma_{3}=\left(\begin{array}[]{cc}0&-i\\ i&0\end{array}\right).$$

$${\mathcal{L}}_{3+1}=\sum _{\sigma }{\overline{\Psi }}_{\sigma }{\Gamma }_{\mu }{D}_{\mu }{\Psi }_{\sigma }+\frac{1}{4}{F}_{\mu \nu }^2+g\mathbf{\varphi }\cdot \overline{\Psi }\frac{\mathbf{\sigma }}{2}\Psi +{\mathcal{L}}_{\mathbf{\varphi }},$$

$$h(\pm \mathbf{Q}+\mathbf{p})=v_F({\mu }_1{p}_x+{\mu }_2{p}_y)+\mathcal{O}({\mathbf{p}}^2),$$

$${\mathcal{L}}_{2+1}$$

$$+{\mathcal{L}}_{{\mathbf{\varphi }}_i},$$

$${\mathcal{L}}_{{\mathbf{\varphi }}_i}$$

$$+\sum _{i<j}\left[u^{\prime }{\mathbf{\varphi }}_i^2{\mathbf{\varphi }}_j^2+u^{\prime \prime }({\mathbf{\varphi }}_i\cdot {\mathbf{\varphi }}_j{)}^2\right]+\mathcal{O}\left(({\mathbf{\varphi }}_i^2{)}^3\right),$$

$${\mathcal{L}}_{2+1}=\sum _{\alpha ,\sigma }{\overline{\psi }}_{\alpha \sigma }{\gamma }_{\mu }{D}_{\mu }{\psi }_{\alpha \sigma }+\frac{1}{4}{F}_{\mu \nu }^2+g\mathbf{\varphi }\cdot {\Psi }^{†}{\gamma }_0\mathbf{\sigma }\Psi +{\mathcal{L}}_{\mathbf{\varphi }},$$

$$\displaystyle\left(\begin{array}[]{c}\psi_{R\sigma}\\ \psi_{L\sigma}\end{array}\right)\rightarrow\left(\begin{array}[]{cc}e^{i\alpha}&0\\ 0&e^{-i\alpha}\end{array}\right)\left(\begin{array}[]{c}\psi_{R\sigma}\\ \psi_{L\sigma}\end{array}\right),$$

$$\mathcal{L}$$

$$+\frac{1}{2}({\partial }_{\mu }\mathbf{\varphi }{)}^2+\frac{1}{2}{m}^2{\mathbf{\varphi }}^2+{\lambda }^2({\mathbf{\varphi }}^2{)}^2+g\mathbf{\varphi }\cdot \sum _{i=1}^{N_f}{\overline{\Psi }}_i\frac{\mathbf{\sigma }}{2}{\Psi }_i,$$

$$\displaystyle\Psi_{i}\equiv\left(\begin{array}[]{c}\Psi_{i\uparrow}\\ \Psi_{i\downarrow}\end{array}\right),$$

$${\mathcal{L}}_0$$

$$+\frac{1}{2}{\left({\partial }_{\mu }{\mathbf{\varphi }}_0\right)}^2+\frac{1}{2}{m}_0^2{\mathbf{\varphi }}_0^2+{\lambda }_0^2{\left({\mathbf{\varphi }}_0^2\right)}^2$$

$$+g_0{\mathbf{\varphi }}_0\cdot \sum _{i=1}^{N_f}{\overline{\Psi }}_i^0\frac{\mathbf{\sigma }}{2}{\Psi }_i^0,$$

$${\mathcal{L}}_R$$

$$+\frac{1}{2}{Z}_{\varphi }({\partial }_{\mu }\mathbf{\varphi }{)}^2+\frac{1}{2}{Z}_{{\varphi }^2}{m}^2{\mu }^2{\mathbf{\varphi }}^2+Z_{{\lambda }^2}{\lambda }^2{\mu }^{ϵ}({\mathbf{\varphi }}^2{)}^2$$

$$+Z_{g}g{\mu }^{ϵ/2}\mathbf{\varphi }\cdot \sum _{i=1}^{N_f}{\overline{\Psi }}_i\frac{\mathbf{\sigma }}{2}{\Psi }_i,$$

$$e^2$$

Full text

I Introduction

The discovery of fractionalized phases of correlated quantum matter which fall outside the standard broken-symmetry classification pioneered by Landau is a prime goal of modern condensed matter physics. The inability to describe such phases in terms of conventional local order parameters suggests that continuous transitions among them, or between them and conventional phases, are themselves exotic and may fall outside the traditional Landau-Ginzburg-Wilson (LGW) paradigm according to which critical properties are solely determined by the long-wavelength, low-energy dynamics of critical order parameter fluctuations.

The paradigmatic example of fractionalized phase of matter is the algebraic spin liquid (ASL) or Dirac spin liquid, originally proposed by Affleck and Marston Affleck and Marston (1988); Marston and Affleck (1989) as a possible ground state of the spin-1/2 Heisenberg model on the 2D square lattice. This state can be envisaged as the result of quantum disordering a magnetically ordered ground state, e.g., the antiferromagnetic (AF) Néel state, in such a way that the bosonic spin-1 magnons of the ordered state fractionalize into a pair of electrically neutral spin-1/2 fermions governed by a Dirac dispersion and interacting with an emergent gauge field. However, it was shown early on Anderson (1952); Huse (1988); Reger and Young (1988); Huse and Elser (1988) that due to the bipartite nature of the square lattice, the nearest-neighbor spin-1/2 Heisenberg model on this lattice develops Néel AF order at zero temperature, and the search for quantum spin liquid ground states rapidly switched its focus to frustrated magnetism Balents (2010), with much emphasis on spin models defined on non-bipartite lattices or with anisotropic interactions.

Surprisingly, a recent sign-problem-free quantum Monte Carlo (QMC) study Xu et al. (2019) demonstrated the realization of an ASL-like phase on the 2D square lattice. The model simulated is a model of spinful fermions coupled to a lattice $U(1)$ gauge field (see Sec. II.1); a gauge constraint (Gauss’ law) enforced dynamically at low temperatures ensures there are no gauge-invariant local fermionic excitations in the low-energy Hilbert space, thus the ground state of such a Hamiltonian corresponds to that of a bosonic spin model. A magnetic flux of $\pi$ per plaquette is spontaneously generated at the mean-field level, similarly to the Affleck-Marston $\pi$-flux phase Flu . A parameter $J$ in the model controls the strength of gauge fluctuations such that for $J$ below a finite critical value there is no symmetry breaking, and AF and valence-bond-solid (VBS) correlations exhibit a power-law decay with distance, as expected for the ASL Kim and Lee (1999); Rantner and Wen (2002); Hermele et al. (2005); hermele2007. The low-energy effective theory of the model in this regime—(2+1)D quantum electrodynamics (QED3) with $N_{f}=1$ flavor of spinful Dirac fermions, where $N_{f}$ is defined precisely here as the number of spin $SU(2)$ doublets of four-component Dirac spinors—is identical to that of the ASL Flu , and for simplicity we will refer to the phase observed numerically as the ASL. At low energies the ASL is believed 111To formally derive the low-energy effective field theory of the ASL one should successively integrate out short-wavelength degrees of freedom from a microscopic spin model, or a lattice gauge theory such as Eq. (II.1), to arrive at a coarse-grained continuum description. However, rigorously implementing such a sequence of exact renormalization group transformations in an interacting many-body system is practically impossible. Generally accepted arguments for the QED3* description starting from a microscopic spin model are presented in Ref. Hermele et al. (2005). Furthermore, numerical studies using variational QMC Iqbal et al. (2011, 2013, 2014) and the density-matrix renormalization group (DMRG) He et al. (2017); Zhu et al. (2018) support the description of the ASL phase in terms of Dirac fermions coupled to a $U(1)$ gauge field to be described by a strongly coupled (2+1)D conformal field theory whose universal properties, notably the critical exponents controlling the power-law correlations mentioned above, can be systematically computed using field-theoretic approaches such as the $\epsilon$-expansion Di Pietro et al. (2016); Di Pietro and Stamou (2017, 2018); Zerf et al. (2018) and the large-$N_{f}$ expansion Gracey (1993a, 1994); Rantner and Wen (2002); Hermele et al. (2005); hermele2007; Chester and Pufu (2016a). As alluded to in the abstract, the ASL phase is characterized by an absence of quasiparticle excitations. As pointed out in the classic papers on the ASL Wen and Lee (1996); Rantner and Wen (2001); Wen (2002); Hermele et al. (2005); hermele2007, the standard low-energy description of the ASL is in terms of Dirac fermions coupled to a fluctuating gauge field, i.e., QED3. At high energies in this theory the coupling is weak, the physical excitations of the system are in one-to-one correspondence with free Dirac fermions, and a quasiparticle description holds. By contrast, at low energies—where the QED3 description of the ASL is believed to hold—the coupling is strong and the quasiparticle description in terms of free Dirac fermions breaks down. Mathematically, spectral functions of physical (i.e., gauge-invariant) observables do not exhibit quasiparticle peaks at low energies, but rather are characterized by power laws with universal exponents.

Remarkably, Ref. Xu et al. (2019) reported not only the first unbiased numerical observation of the ASL, but also that of a direct quantum phase transition from the ASL to the AF Néel state at a critical value of the parameter $J$ introduced above. (The existence of this transition is nontrivial because no intermediate phases appear, by contrast with the case with twice as many fermion flavors where the ASL-like phase first transitions to a VBS phase, then to an AF phase.) Numerical evidence further points to a continuous transition, as opposed to a first-order transition, implying the existence of universal critical phenomena. On general grounds one expects the low-energy theory of the transition to be of the QED3-Gross-Neveu-Yukawa (GNY) type Lu et al. (2017); Xu et al. (2019); Song et al. (2019), as for the recently studied transitions between the ASL and a gapped chiral spin liquid Janssen and He (2017); Ihrig et al. (2018a); Zerf et al. (2018); Boyack et al. (2019) or a gapped $\mathbb{Z}_{2}$ spin liquid Boyack et al. (2018), but with the spontaneous breaking of an $O(3)$ symmetry as appropriate for Néel order Ghaemi and Senthil (2006), as opposed to the Ising and XY symmetries appropriate for those two transitions, respectively. In theories of this type, gapless $U(1)$ gauge fluctuations as well as gapless Dirac fermions couple strongly to bosonic order parameter fluctuations and can give rise to novel universality classes of critical behavior distinct from the standard (Wilson-Fisher) Ising, XY, and Heisenberg universality classes, provided bona fide critical fixed points of the renormalization group (RG) exist in those theories. While a stable critical point was found for the chiral Ising QED3-GNY model both in the $\epsilon$-expansion Janssen and He (2017); Ihrig et al. (2018a); Zerf et al. (2018) for all $N_{f}$ and in the large-$N_{f}$ expansion Gracey (1992, 1993b); *gracey2018b; Gracey (1993c, 2018b); Boyack et al. (2019), in the chiral XY QED3-GNY model Boyack et al. (2018) a fixed point found at one-loop order in the $\epsilon$ expansion was shown to disappear below a certain critical value of $N_{f}$, in analogy to the phenomenon of fluctuation-induced first-order transitions Halperin et al. (1974), while the large-$N_{f}$ expansion predicts a stable critical point. Whether the chiral Heisenberg QED3-GNY model supports a critical point that remains stable at all $N_{f}$, and if so, what its universal critical properties are, are thus important open questions in light of the numerical discoveries reported in Ref. Xu et al. (2019). In Ref. Ghaemi and Senthil (2006), the chiral Heisenberg QED3-GNY model was proposed as the critical theory for a transition between the algebraic spin liquid, in its staggered flux phase realization on the square lattice, and a Néel antiferromagnet. The authors studied this model in the $\epsilon$ expansion below four dimensions at one-loop order, for the specific case of $N_{f}=1$, and found a critical fixed point describing a continuous phase transition. However, as we will show below, their RG beta functions and critical exponents disagree with ours. Our combined $\epsilon$-expansion and large-$N_{f}$ analyses allows us to perform several nontrivial consistency checks, detailed below, which lead us to believe our results are correct.

In this work, we revisit the problem of the Néel–algebraic-spin-liquid transition, motivated by the numerical results of Ref. Xu et al. (2019), and present a detailed study of the critical properties of the chiral Heisenberg QED3-GNY model for a generic number $N_{f}$ of flavors of four-component Dirac fermions carrying an additional $SU(2)$ spin index. After deriving the specific form of the low-energy Lagrangian from microscopic considerations (Sec. II), we employ two complementary approaches to study the critical regime: an $\epsilon$ expansion in $d=4-\epsilon$ spacetime dimensions carried out to four-loop order for arbitrary $N_{f}$ (Sec. III), as well as a large-$N_{f}$ expansion carried out up to second order in fixed but arbitrary $d$ dimensions (Sec. IV). The use of arbitrary $d$ in the latter approach allows us to establish consistency between the two methods, and gives us confidence in the correctness of our analysis. We do, however, pay attention to contributions specific to $d=3$ (Aslamazov-Larkin diagrams) which arise for certain quantities in the large-$N_{f}$ expansion. Through the $\epsilon$ expansion we establish the existence of a stable critical fixed point for all $N_{f}$, including the $N_{f}=1$ case relevant for the Néel-ASL transition on the 2D square lattice observed in the QMC study Xu et al. (2019). Gauge-invariant universal critical properties including the order parameter anomalous dimension $\eta_{\phi}$, the (inverse) correlation length exponent $1/\nu$, and the stability critical exponent $\omega$, but also the universal exponents $\Delta_{\text{CDW}}$, $\Delta_{\text{VBS}}$, and $\Delta_{\text{QAH}}$ characterizing the power-law decay of charge-density-wave (CDW), VBS, and quantum anomalous Hall (QAH) correlations at the transition, are computed to $\mathcal{O}(\epsilon^{4})$ and up to $\mathcal{O}(1/N_{f}^{2})$. While the CDW, VBS, and QAH exponents have been computed in the ASL phase before Hermele et al. (2005); *hermele2007, the Néel-ASL critical point and its arbitrary-$N_{f}$ generalization correspond to a different (2+1)D conformal field theory than that describing the ASL phase, and is thus characterized by different—but equally universal—values for these exponents. We then use Padé and Padé-Borel resummation techniques to obtain numerical estimates of the critical exponents (Sec. V). As a byproduct of our calculations, we also obtain CDW and VBS exponents at the semimetal-AF insulator quantum phase transition seen in the Hubbard model on the honeycomb lattice Sorella and Tosatti (1992); Herbut (2006); Honerkamp (2008); Herbut et al. (2009); Sorella et al. (2012); Ulybyshev et al. (2013); Assaad and Herbut (2013); Janssen and Herbut (2014); Parisen Toldin et al. (2015); Otsuka et al. (2016); Buividovich et al. (2018) (see Sec. III.7 and IV.4). These correspond to the scaling dimensions of certain fermion bilinear operators in the chiral Heisenberg GNY model, which to our knowledge have not been computed before despite recent studies of this model Zerf et al. (2017); Gracey (2018c).

II The Néel–algebraic-spin-liquid transition

In this section we derive the continuum quantum field theory describing the Néel-ASL transition on the square lattice (Sec. II.1) as observed in the QMC study of Ref. Xu et al. (2019), and contrast it with the field theory description of similar magnetic ordering transitions on the kagomé lattice (Sec. II.2).

II.1 Square lattice

The starting point is the compact $U(1)$ lattice gauge theory model considered in Ref. Xu et al. (2019), here restricted to two fermion flavors:

[TABLE]

Here $c_{r\sigma}^{({\dagger})}$ is the annihilation (creation) operator for a fermion with spin $\sigma=\uparrow,\downarrow$ on site $r$ of a 2D square lattice, $\theta_{r,r+\hat{\mu}}\in[0,2\pi)$ is an angular variable living on the nearest-neighbor links of the square lattice, with $\mu=x,y$ denoting its two orthogonal directions, and $L_{r,r+\hat{\mu}}$ is the operator conjugate to $\theta_{r,r+\hat{\mu}}$, i.e., $[\theta_{r,r+\hat{\mu}},L_{r^{\prime},r^{\prime}+\hat{\nu}}]=i\delta_{rr^{\prime}}\delta_{\mu\nu}$. Both $\theta_{r,r+\hat{\mu}}$ and $L_{r,r+\hat{\mu}}$ can be interpreted as vector fields $\theta_{r,\mu}$ and $L_{r,\mu}$, and $\Delta\times\bm{\theta}=\Delta_{x}\theta_{r,y}-\Delta_{y}\theta_{r,x}$ denotes the lattice curl of $\theta_{r,\mu}$, corresponding to its circulation around an elementary plaquette $\square$, where $\Delta_{\mu}\theta_{r,\nu}=\theta_{r+\hat{\mu},\nu}-\theta_{r,\nu}$ is the lattice derivative. The first term, proportional to $J$, is an electric-field term that controls the strength of the gauge fluctuations, while the last term, proportional to $K>0$, is a magnetic-field term that favors the maximal (i.e., $\pi$) amount of magnetic flux per plaquette. The fermion chemical potential is set to zero, corresponding to half filling.

When $J=0$, $\theta_{r,\mu}$ is a classical variable and the fermions are noninteracting. The ground state at half-filling is guaranteed by Lieb’s theorem Lieb (1994) to have $\pi$ flux per plaquette. When $J$ is nonzero but small, the gauge field acquires dynamics but the QMC results Xu et al. (2019) indicate the $\pi$-flux phase remains stable up to some critical value, given by $J_{c}=1.43(46)$ for $t=K=1$. In the presence of gauge fluctuations the $\pi$-flux phase is an ASL described by the conformally invariant, infrared-stable fixed point of QED3 with two flavors of four-component fermions. For $J>J_{c}$ the model develops long-range $(\pi,\pi)$ AF order, consistent with the known Néel ground state of the $SU(2)$-symmetric Heisenberg model, to which model (II.1) reduces in the $J\rightarrow\infty$ limit. The smooth behavior at $J=J_{c}$ of both the AF correlation ratio and the flux energy per plaquette observed in the numerical results Xu et al. (2019) suggests the Néel-ASL transition is continuous, and should be described by a critical fixed point of a continuum quantum field theory. To derive the precise form of this quantum field theory, which is expected to be of the chiral Heisenberg QED3-GNY type, we simply need to express the Néel order parameter in terms of the low-energy Dirac fermions in the $\pi$-flux phase. This in turn will dictate the correct form of the Yukawa coupling between Dirac fermions and the $O(3)$ vector field describing order parameter fluctuations. As will be seen in detail below and in Sec. II.2, owing to the precise relation between the low-energy Dirac fermion fields and the original microscopic lattice fermions, which differs on the square and kagomé lattices, the (2+1)D Néel-ASL transition on the square lattice can be accessed by dimensional continuation from the four-dimensional Lorentz-invariant chiral Heisenberg QED3-GNY model, while conceptually similar magnetic ordering transitions on the kagomé lattice cannot. This has important implications for $\epsilon$-expansion studies of such transitions as will be discussed below. Similar difficulties in dimensional continuation for critical points described by the chiral Ising QED3-GNY model were discussed in Ref. Boyack et al. (2019), but in a purely field-theoretical context without reference to a microscopic lattice. Related issues were discussed recently for pure QED3 in the high-energy literature Di Pietro et al. (2016); Di Pietro and Stamou (2017), and for the pure chiral XY GNY model in Ref. Wamer and Affleck (2018).

We begin by studying the $\pi$-flux phase in the absence of gauge fluctuations, described by the fermionic part of the $J=0$ limit of the model (II.1). Choosing a gauge such that the hopping amplitude on alternating links of the square lattice acquires an additional phase factor of $-1$ (Fig. 1), the Hamiltonian can be written as

[TABLE]

where $c_{{\bm{k}}\sigma}=(c_{{\bm{k}},A,\sigma},c_{{\bm{k}},B,\sigma})$ is a two-component spinor, $A,B$ denote the two sublattices of the square lattice, and

[TABLE]

where

[TABLE]

The integral is over wavevectors ${\bm{k}}$ (measured in inverse units of the lattice constant) in the first Brillouin zone, which is reduced due to the doubling of the unit cell stemming from the choice of gauge and forms a square of side length $\pi\sqrt{2}$ rotated by 45∘ with respect to the original Brillouin zone. The spectrum is given by $E_{\pm}({\bm{k}})=\pm 2t\sqrt{\sin^{2}k_{x}+\cos^{2}k_{y}}$ and exhibits two inequivalent Dirac points at $\pm{\bm{Q}}$ with ${\bm{Q}}\equiv\left(0,\frac{\pi}{2}\right)$. At half filling the low-energy excitations are particle-hole excitations of the Dirac fermions, and we can linearize the Hamiltonian near the Dirac points,

[TABLE]

where $|{\bm{p}}|$ is much less than the dimensions of the Brillouin zone, we denote Pauli matrices acting on the Dirac components by $\mu_{1},\mu_{2},\mu_{3}$, and $v_{F}\equiv 2t$. The low-energy Hamiltonian is thus given by

[TABLE]

where we define the slow (Dirac) fields

[TABLE]

Note that defining $\psi_{-\sigma}({\bm{p}})$ as $i\mu_{2}$ acting on the two-component spinor $c_{-{\bm{Q}}+{\bm{p}},\sigma}$ is necessary to absorb the extra minus sign in front of $\mu_{1}$ in the Hamiltonian matrix (7) linearized near $-{\bm{Q}}$. Introducing the two-component Dirac adjoint $\overline{\psi}_{\alpha\sigma}\equiv\psi^{\dagger}_{\alpha\sigma}\gamma_{0}$ with $\gamma_{0}=\mu_{3}$, and Fourier transforming $H_{0}$ back to (continuous) real space, we find that the imaginary-time Lagrangian density $\mathcal{L}_{0}=\sum_{\alpha,\sigma}\psi_{\alpha\sigma}^{\dagger}\partial_{0}\psi_{\alpha\sigma}+H_{0}$ can be written as

[TABLE]

where $(\gamma_{0},\gamma_{1},\gamma_{2})=(\mu_{3},\mu_{2},-\mu_{1})$ are $2\times 2$ Euclidean gamma matrices, and we have rescaled the spatial coordinates to set the Dirac velocity $v_{F}$ to unity.

We now turn to Néel order. Referring to Fig. 1, the Néel order parameter ${\bm{N}}$ is an $SU(2)$ triplet of fermion bilinears given by

[TABLE]

where $\mathcal{N}$ is the total number of sites, ${\bm{R}}=n_{1}{\bm{a}}_{1}+n_{2}{\bm{a}}_{2}$ with $n_{1},n_{2}\in\mathbb{Z}$ denotes sites of the underlying Bravais lattice with primitive lattice vectors ${\bm{a}}_{1}\equiv\hat{{\bm{x}}}-\hat{{\bm{y}}}$ and ${\bm{a}}_{2}\equiv\hat{{\bm{x}}}+\hat{{\bm{y}}}$, measured in units of the lattice constant, and ${\bm{S}}_{{\bm{R}}A(B)}=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}c_{{\bm{R}}A(B)\sigma}^{\dagger}\bm{\sigma}_{\sigma\sigma^{\prime}}c_{{\bm{R}}A(B)\sigma^{\prime}}^{\phantom{\dagger}}$ is the local spin operator on sublattice $A(B)$, where $\bm{\sigma}=(\sigma_{1},\sigma_{2},\sigma_{3})$ denotes Pauli matrices acting on the (physical) spin components. Fourier transforming and keeping only the low-energy Dirac degrees of freedom near $\pm{\bm{Q}}$, we find

[TABLE]

where we have defined the four-component Dirac spinors

[TABLE]

This motivates the following choice of (reducible) four-dimensional representation of the $SO(3)$ Clifford algebra,

[TABLE]

and the corresponding definition of four-component Dirac adjoint $\overline{\Psi}_{\sigma}\equiv\Psi_{\sigma}^{\dagger}\Gamma_{0}$, in terms of which the Néel order parameter is written as

[TABLE]

with trace over spin indices understood.

Néel order is not explicit in the microscopic Hamiltonian (II.1) but arises from the spontaneous breaking of its spin rotation symmetry. Having identified the form of the Néel order parameter bilinear in terms of the low-energy Dirac fields, we can write down an effective LGW Lagrangian valid near the Néel-ASL transition in physical 2+1 dimensions,

[TABLE]

where $\bm{\phi}$ is a Lorentz-scalar $O(3)$ vector field that has the interpretation of dynamical order parameter field, and $\mathcal{L}_{\bm{\phi}}$ is the Lagrangian controlling its dynamics, i.e., the $O(3)$ vector model. In the language of this effective field theory, the ASL phase corresponds to a phase in which $\langle{\bm{\phi}}\rangle=0$ and the $O(3)$ symmetry is unbroken, while in the Néel phase $\langle{\bm{\phi}}\rangle\neq 0$ and the $O(3)$ symmetry is spontaneously broken to an $O(2)$ subgroup of spin rotations about the direction of $\langle{\bm{\phi}}\rangle$. We have accounted for the general situation with gauge fluctuations ($J\neq 0$) by replacing $\partial_{\mu}$ in Eq. (14) by the gauge-covariant derivative $D_{\mu}=\partial_{\mu}-ieA_{\mu}$, with $e$ the gauge coupling and $A_{\mu}$ the continuum limit of the $U(1)$ lattice gauge field $\theta_{r,r+\hat{\mu}}$. We have also added a Maxwell term for the gauge field with $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ the field-strength tensor, which originates microscopically from the terms proportional to $J$ and $K$ in the lattice Hamiltonian (II.1). The Lagrangian (26) is equivalent to that studied in Ref. Ghaemi and Senthil (2006), although the fermionic mean-field state considered in the latter reference is the staggered flux state rather than the $\pi$-flux state, and thus the projective symmetry group characterizing the corresponding spin liquid state is different (see Appendix A).

As in previous work on related theories Di Pietro et al. (2016); Di Pietro and Stamou (2017, 2018); Janssen and He (2017); Ihrig et al. (2018a); Zerf et al. (2018); Boyack et al. (2018), we will be interested in studying the critical properties of Eq. (26) in the standard $\epsilon$ expansion below four (spacetime) dimensions, and thus require a four-dimensional Lorentz-invariant Lagrangian which dimensionally continues to Eq. (26). In $d=4$ dimensions the reducible four-dimensional representation of the $SO(3)$ Clifford algebra becomes an irreducible representation of the $SO(4)$ Clifford algebra when supplemented with a fourth gamma matrix $\Gamma_{3}$, which for the choice of $SO(3)$ gamma matrices in Eq. (24) is given by

[TABLE]

One can then check explicitly that the four-dimensional Lagrangian

[TABLE]

is invariant under four-dimensional Lorentz transformations, and reduces to the physical Lagrangian (26) when all terms involving $x_{3}$ and $A_{3}$ are ignored. Thus the critical properties of the Néel-ASL transition on the square lattice can be consistently studied in an $\epsilon$ expansion of the chiral Heisenberg QED-GNY theory (30) with a single $SU(2)$ doublet $(\Psi_{\uparrow},\Psi_{\downarrow})$ of four-component Dirac fermions.

Finally, we note that while our field-theoretic study is motivated by the numerical study of a specific model, Hamiltonian (II.1), it is more broadly applicable as it relies only on the universal low-energy description of the ASL phase in terms of QED3, and the broken symmetries in the Néel phase. In other words, if a Néel-ASL transition is found in any other microscopic spin model on the square lattice, its universal critical exponents would be given by those we have calculated, regardless of the particular microscopic model. We choose to focus on the model (II.1) as it is the only microscopic model in which a direct, continuous transition from the ASL to the Néel state has been observed in unbiased, controlled QMC simulations Xu et al. (2019).

II.2 Kagomé lattice

We now briefly review the slave-fermion description of the ASL on the kagomé lattice Hastings (2000); Ran et al. (2007); Hermele et al. (2008) and its possible magnetic ordering instabilities. Similar to the $\pi$-flux phase on the square lattice, ignoring gauge fluctuations the ASL on the kagomé lattice can be understood as the ground state of noninteracting fermions with nearest-neighbor hopping on the kagomé lattice and $\pi$ flux through each hexagonal plaquette, but zero flux through each triangular plaquette. The simplest gauge choice doubles the unit cell in (say) the ${\bm{a}}_{1}$ direction, where ${\bm{a}}_{1}=(2,0)$ and ${\bm{a}}_{2}=(1,\sqrt{3})$ are the primitive vectors of the underlying triangular Bravais lattice, with lengths measured in units of the side length of a kagomé triangle. Correspondingly the original (hexagonal) first Brillouin zone is halved (and becomes rectangular); at half filling the Fermi energy intersects two inequivalent Dirac points at $\pm{\bm{Q}}=\frac{\pi}{2\sqrt{3}}(0,\pm 1)$. Effective Hamiltonians similar to Eq. (7) can be obtained in first-order degenerate perturbation theory,

[TABLE]

with $v_{F}=\sqrt{2}t$. Slow Dirac fields can be defined in analogy to Eq. (13), albeit with a more complicated dependence on the microscopic lattice fields (the doubled unit cell has six sites), and the imaginary-time Lagrangian density is again given by Eq. (14) with the same set of $2\times 2$ gamma matrices, $(\gamma_{0},\gamma_{1},\gamma_{2})=(\mu_{3},\mu_{2},-\mu_{1})$ Hermele et al. (2008). As before, gauge fluctuations are incorporated by promoting $\partial_{\mu}$ to the gauge-covariant derivative $D_{\mu}$ in the Dirac Lagrangian and adding a Maxwell term for the emergent continuum gauge field $A_{\mu}$.

We now turn to magnetic ordering instabilities, i.e., instabilities accompanied by spontaneous breaking of the spin $SU(2)$ symmetry (which excludes VBS states) of an underlying microscopic spin model on the kagomé lattice, such as the spin-1/2 Heisenberg antiferromagnet Hastings (2000); Ran et al. (2007); Hermele et al. (2008); Iqbal et al. (2011, 2013, 2014); He et al. (2017); Zhu et al. (2018). We further focus on patterns of ordering that can be achieved through the condensation of Dirac fermion bilinears which gaps out the Dirac fermions (once the fermions are gapped, monopole proliferation typically ensues, which further gaps out the dynamical photon and leads to confinement Song et al. (2019, 2018)—see Sec. IV.5). This leaves two possibilities Hermele et al. (2008): a three-fold degenerate collinear stripe AF state with ordering wavevectors corresponding to the three symmetry-related $M$ points of the physical (gauge-invariant) hexagonal Brillouin zone, and a zero-wavevector noncollinear state with nonzero vector spin chirality Lu et al. (2017).

We first consider the stripe AF state (i.e., Fig. 4 in Ref. Hermele et al. (2008)), which corresponds to the condensation of either of the three degenerate spin-triplet, time-reversal-odd bilinears $\Psi^{\dagger}\mu_{3}\tau_{i}\bm{\sigma}\Psi$, $i=1,2,3$, where $\tau_{1},\tau_{2},\tau_{3}$ are Pauli matrices acting on the nodal/valley index $\alpha=\pm$ distinguishing the two Dirac nodes Hermele et al. (2008). Introducing a corresponding set of Lorentz-scalar $O(3)$ vector fields $\bm{\phi}_{1},\bm{\phi}_{2},\bm{\phi}_{3}$ describing the fluctuating order parameters, the form of the LGW Lagrangian for the stripe AF-ASL transition is dictated by the transformation properties of the stripe AF fermion bilinears under space group operations, which have been worked out in detail in Ref. Hermele et al. (2008). One obtains

[TABLE]

where the order parameter Lagrangian $\mathcal{L}_{\bm{\phi}_{i}}$ is given by

[TABLE]

where $r,u,u^{\prime},u^{\prime\prime},\ldots$ are phenomenological coupling constants. The $\mathbb{Z}_{3}$ anisotropy introduced by the quartic terms in Eq. (II.2) is a consequence of the discrete $C_{6}$ symmetry of the kagomé lattice which rotates the three $M$ points into each other.

In order to study the critical properties of the stripe AF-ASL transition in the $4-\epsilon$ expansion one must lift the (2+1)D theory (II.2) to a four-dimensional Lorentz-invariant theory with four-component Dirac spinors. The main issue, as compared to the theory (26) for the Néel-ASL transition on the square lattice, is that the Yukawa coupling involves the nodal matrices $\tau_{i}$. These matrices generate an $SU(2)_{\text{nodal}}$ subgroup of the enlarged $SU(4)$ flavor symmetry of the pure ASL state, i.e., the internal symmetry of the first term in Eq. (II.2), but in the four-dimensional theory elements of this subgroup act as Lorentz group elements on four-component Dirac spinors. In other words, the Yukawa coupling breaks the $SU(2)_{\text{nodal}}$ symmetry in the (2+1)D theory, which means it breaks Lorentz symmetry in the four-dimensional theory. For instance, suppose we choose the $4\times 4$ gamma matrix representation given in Eq. (24) and Eq. (29): then the three bilinears $\Psi^{\dagger}\gamma_{0}\tau_{i}\bm{\sigma}\Psi$, $i=1,2,3$ are proportional to $i\overline{\Psi}\Gamma_{3}\bm{\sigma}\Psi$, $-i\overline{\Psi}\Gamma_{5}\bm{\sigma}\Psi$, and $\overline{\Psi}\bm{\sigma}\Psi$, respectively, only the last two of which are invariant under proper orthochronous Lorentz transformations (recall that $\Gamma_{5}\equiv\Gamma_{0}\Gamma_{1}\Gamma_{2}\Gamma_{3}$ is a Lorentz pseudoscalar). Different choices of $4\times 4$ gamma matrices that preserve the equality between the first term in Eq. (II.2) and its four-component counterpart (i.e., the first term of Eq. (30)) in the three-dimensional limit simply permute the assignments of $\Psi^{\dagger}\gamma_{0}\tau_{i}\bm{\sigma}\Psi$ to those same three bilinears.

A slightly more general argument of the above conclusion is as follows. There are three independent stripe AF ground states related by point-group rotations, described in the field theory by the condensation of three fermion bilinears. For such bilinears to open a gap in the Dirac fermion spectrum, they must be mass terms, i.e., they must be given by linearly independent matrices which anticommute with the massless Dirac Hamiltonian. In 3+1 dimensions, there are only two such Lorentz-invariant mass terms: the normal mass $\overline{\Psi}\Psi$ and the axial mass $i\overline{\Psi}\Gamma_{5}\Psi$. (Since AF order breaks spin rotation symmetry, the appropriate mass terms are actually $\overline{\Psi}\bm{\sigma}\Psi$ and $i\overline{\Psi}\Gamma_{5}\bm{\sigma}\Psi$, but spin is an internal degree of freedom which does not play a role in this argument.) Therefore, one of the three mass terms describing stripe AF order must necessarily violate Lorentz invariance.

We finally turn to the ${\bm{q}}=0$ noncollinear state with vector spin chirality order (i.e., Fig. 7 in Ref. Hermele et al. (2008)), which corresponds to the condensation of a unique spin-triplet bilinear $\Psi^{\dagger}\mu_{3}\bm{\sigma}\Psi$ even under time reversal Hermele et al. (2008). The LGW theory for the transition from the ASL to this state was given in Ref. Lu et al. (2017),

[TABLE]

and involves a single $O(3)$ vector $\bm{\phi}$. In four-component notation, the bilinear appearing in the Yukawa coupling is proportional to $i\overline{\Psi}\Gamma_{3}\Gamma_{5}\bm{\sigma}\Psi$ for standard choices of $4\times 4$ gamma matrices that preserve the form of the pure ASL Lagrangian; in such cases the theory (34) thus breaks Lorentz invariance when lifted to four dimensions.

To demonstrate generally that the lift from $d=2+1$ to $d=3+1$ for the ${\bm{q}}=0$ noncollinear order requires breaking Lorentz invariance is more subtle, as there is a unique ground state (modulo global $SU(2)$ spin rotations), which implies only one mass term is required. In general, one must consider the fate of the theory’s symmetries upon dimensional continuation. In $d=2+1$, the ASL has an emergent $SU(4)\supset SU(2)_{\text{nodal}}\times SU(2)_{\text{spin}}$ “flavor” symmetry which rotates the four two-component Dirac spinors $(\psi_{R\uparrow},\psi_{R\downarrow},\psi_{L\uparrow},\psi_{L\downarrow})$ into one another. (Here we denote the nodal indices by $R,L$ instead of $+,-$ used previously, for reasons that will become clear shortly.) The spin subgroup $SU(2)_{\text{spin}}$ acts only on the spin indices $\uparrow,\downarrow$ while the nodal subgroup $SU(2)_{\text{nodal}}$ acts only on the nodal indices $R,L$. In QED4, the two-component Dirac fermions become two-component Weyl fermions: $\psi_{R\sigma}$ forms an $SU(2)$ doublet of two-component right-handed Weyl fermions, while $\psi_{L\sigma}$ forms an $SU(2)$ doublet of left-handed Weyl fermions. (The two Weyl components must have opposite chirality, otherwise the $d=3+1$ theory would have a gauge anomaly, which would make it inconsistent.) Among other symmetries, massless QED4 has a $U(1)_{R-L}$ chiral symmetry, under which the Weyl spinors transform as $\psi_{R\sigma}\rightarrow e^{i\alpha}\psi_{R\sigma}$ and $\psi_{L\sigma}\rightarrow e^{-i\alpha}\psi_{L\sigma}$. It is well known that this chiral symmetry protects the masslessness of the four-component Dirac fermion $\Psi_{\sigma}=(\psi_{R\sigma},\psi_{L\sigma})$. (In a condensed matter setting, this is simply the topological protection of Weyl semimetals.) Conversely, mass terms—such as the only two Lorentz-invariant mass terms, $\overline{\Psi}\Psi$ and $i\overline{\Psi}\Gamma_{5}\Psi$—couple the right-handed and left-handed Weyl components and break this chiral symmetry. Now, the $U(1)_{R-L}$ transformation can be written as

[TABLE]

which is an $SU(2)_{\text{nodal}}$ rotation in the (2+1)-dimensional interpretation. However, the bilinear $\Psi^{\dagger}\mu_{3}\bm{\sigma}\Psi$ corresponding to the ${\bm{q}}=0$ noncollinear order is proportional to the identity in nodal space (i.e., it does not contain nodal Pauli matrices $\tau_{i}$); by contrast with Néel order on the square lattice, it produces the same mass (amplitude and sign) at both Dirac nodes. In other words, this bilinear preserves the $SU(2)_{\text{nodal}}\supset U(1)_{R-L}$ symmetry, and thus cannot correspond to a Lorentz-invariant mass term in 3+1 dimensions.

III $\epsilon$-expansion

We have shown above that the LGW theories for transitions from the ASL to either the stripe AF or the ${\bm{q}}=0$ noncollinear state on the kagomé lattice break Lorentz invariance when lifted from three to four spacetime dimensions. This creates difficulties when applying the standard $\epsilon$ expansion as a variety of Lorentz-breaking terms not present in the original Lagrangian will be generated under renormalization (for a recent discussion of similar issues in the pure GNY model, see Ref. Wamer and Affleck (2018)). By contrast, the LGW theory (26) for the Néel-ASL transition on the square lattice can be extended to a four-dimensional Lorentz-invariant theory, Eq. (30). This allows for an RG analysis in the usual $\epsilon$ expansion, to which we now turn.

III.1 Renormalization group

The starting point is the Lagrangian of the four-dimensional chiral Heisenberg QED-GNY model,

[TABLE]

where $\not{D}\equiv\Gamma_{\mu}(\partial_{\mu}-ieA_{\mu})$ is the gauge-covariant derivative and $\Gamma_{\mu}$ are $4\times 4$ Euclidean gamma matrices representing the $SO(4)$ Clifford algebra $\{\Gamma_{\mu},\Gamma_{\nu}\}=2\delta_{\mu\nu}$. We denote by $N_{f}$ the number of flavors of $SU(2)$ doublets of four-component Dirac fermions, i.e.,

[TABLE]

where $\Psi_{i\sigma}$, $\sigma=\uparrow,\downarrow$ are four-component Dirac fermions. With this convention the Lagrangian (30) for the Néel-ASL transition on the square lattice corresponds to $N_{f}=1$. The triplet of real fields $\bm{\phi}=(\phi_{1},\phi_{2},\phi_{3})$ transforms in the vector representation of $SU(2)$, and $\xi$ is a gauge-fixing parameter ($\xi=0$ is the Landau gauge). By contrast, the chiral Ising QED-GNY model studied in Ref. Janssen and He (2017); Ihrig et al. (2018a); Zerf et al. (2018) possesses a single real scalar field coupled to a fermion mass term, instead of a spin bilinear, and the chiral XY QED-GNY model studied in Ref. Boyack et al. (2018) possesses a complex scalar field coupled to a superconducting mass term.

We study the critical (i.e., at $m^{2}=0$) properties of the model (III.1) in $d=4-\epsilon$ spacetime dimensions with field-theoretic RG using the modified minimal subtraction ($\overline{\text{MS}}$) prescription. We consider a bare Lagrangian,

[TABLE]

in terms of bare fields $\Psi_{i}^{0}$, $\bm{\phi}_{0}$, $A_{\mu}^{0}$ and coupling constants $e_{0}$, $\xi_{0}$, $m_{0}$, $\lambda_{0}$, $g_{0}$, with $\not{D}^{0}=\Gamma_{\mu}(\partial_{\mu}-ie_{0}A_{\mu}^{0})$, and a renormalized Lagrangian

[TABLE]

where $\not{D}=\Gamma_{\mu}(\partial_{\mu}-ie\mu^{\epsilon/2}A_{\mu})$. We define renormalized fields $\Psi_{i}=Z_{\Psi}^{-1/2}\Psi_{i}^{0}$, $A_{\mu}=Z_{A}^{-1/2}A_{\mu}^{0}$, $\bm{\phi}=Z_{\phi}^{-1/2}\bm{\phi}_{0}$ and dimensionless renormalized coupling constants,

[TABLE]

where $\mu$ is an arbitrary renormalization scale. The renormalization constants $Z_{X}$, $X=\Psi,A,\phi,\phi^{2},\lambda^{2},g$ are calculated up to four-loop order using an automated setup, technical details of which can be found in previous publications Zerf et al. (2016); Mihaila et al. (2017); Zerf et al. (2017, 2018). The number of diagrams that arise during the perturbative calculation of the renormalization constants is exactly the same as in the chiral-Ising GNY-QED3 case Zerf et al. (2018), because one needs only to “dress” the chiral-Ising Lorentz amplitudes with $SU(2)$ weight-factors.

III.2 Beta functions

The flow of the running couplings $e^{2}$, $g^{2}$, $\lambda^{2}$ is governed by the following RG equations,

[TABLE]

where the beta functions are defined as

[TABLE]

for $\alpha=e,g,\lambda$, and we absorb a factor of $1/(4\pi)^{2}$ in the definition of the couplings, $\alpha^{2}/(4\pi)^{2}\rightarrow\alpha^{2}$. Likewise, the anomalous dimensions $\gamma_{X}$ are defined as

[TABLE]

Note that Eq. (53)-(55) have the same form as for the chiral Ising QED-GNY model Zerf et al. (2018), but the explicit form of the renormalization constants, and thus the anomalous dimensions and beta functions, are different for the two theories. As in Ref. Zerf et al. (2018) we express the beta functions as a sum over contributions at fixed loop order,

[TABLE]

for $\alpha=e,g,\lambda$. We give here the beta functions up to and including three-loop order; contributions at four-loop order are lengthy and made available in electronic format Sup . The beta function for the gauge coupling $e^{2}$ is given by

[TABLE]

The beta function for the Yukawa coupling $g^{2}$ is given by

[TABLE]

Finally, the beta function for the four-scalar coupling $\lambda^{2}$ is given by

[TABLE]

In those expressions $\zeta_{z}$ denotes the Riemann zeta function.

The beta functions above can be checked against previous results in various limits. In the limit $g^{2}=\lambda^{2}=0$ the model reduces to pure QED with $2N_{f}$ flavors of four-component Dirac fermions, and Eq. (59)-(61) together with the four-loop contribution Sup reproduce the four-loop QED beta function Gorishny et al. (1991). For $e^{2}=g^{2}=0$ the model reduces to the bosonic $O(3)$ vector model; Eq. (65)-(67) with the four-loop contribution Sup agree with the known four-loop result vla . Setting $e^{2}=0$ only, one obtains the ungauged chiral Heisenberg GNY model; our result for $\beta_{g^{2}}$ and $\beta_{\lambda^{2}}$ agrees in that limit with the four-loop beta functions recently computed for this model Zerf et al. (2017).

One-loop beta functions for the specific case of $N_{f}=1$, with all couplings $e^{2}$, $g^{2}$, $\lambda^{2}$ nonzero, were obtained in Ref. Ghaemi and Senthil (2006). Modulo trivial rescalings of the coupling constants, our result for $\beta_{e^{2}}$ agrees with theirs, but $\beta_{g^{2}}$ and $\beta_{\lambda^{2}}$ disagree. In Sec. IV, we show that the critical exponents obtained from our beta functions agree order by order up to $\mathcal{O}(\epsilon^{4},1/N_{f}^{2})$ with continuous-$d$ results obtained in the large-$N_{f}$ expansion. We are thus confident that our beta functions, and the critical exponents we derive from them in the following sections, are correct.

III.3 Quantum critical point and stability exponent

In this section we explore the Néel-ASL quantum critical point (QCP), and we will find it has different critical exponents from the ASL phase, which is described by conformal QED3. The ASL phase and the Néel-ASL QCP are described by different conformal field theories; these correspond to different RG fixed points, characterized by a different type of stability and different critical exponents. The ASL phase Hermele et al. (2005) is an example of a critical phase, that is, a gapless phase of matter with power-law correlations and no relevant symmetry-preserving perturbations (i.e., it is a stable phase of matter). A QCP is also characterized by power-law correlations, but is not a stable phase, i.e., it has a relevant symmetry-preserving perturbation corresponding to the tuning parameter for the transition. To summarize, the fact that the exponents are different in the ASL phase and at the QCP is a generic feature of continuous transitions out of critical phases. (As another example, in the context of the AF quantum phase transition for repulsively interacting fermions on the honeycomb lattice, critical exponents go from the free-Dirac universality class in the semimetallic phase to the chiral Heisenberg GNY universality class at the QCP Herbut (2006).) This can be understood physically from the fact that additional soft degrees of freedom—order parameter fluctuations—appear at the QCP, and change the universality class from that of the critical phase itself.

To identify a possible QCP we look for stable fixed points of the RG flow on the critical ($m^{2}=0$) hypersurface, i.e., common zeros of the set of beta functions (53)-(55). We denote fixed-point couplings by $(e^{2}_{*},g^{2}_{*},\lambda^{2}_{*})$. At one-loop order we find eight fixed points: the Gaussian fixed point $(0,0,0)$, the conformal QED fixed point $(\frac{3\epsilon}{16N_{f}},0,0)$, the $O(3)$ Wilson-Fisher fixed point $(0,0,\frac{\epsilon}{88})$, a conformal QED $\times$ $O(3)$ Wilson-Fisher fixed point $(\frac{3\epsilon}{16N_{f}},0,\frac{\epsilon}{88})$, two $O(3)$ GNY fixed points with $e^{2}_{*}=0$ and $g^{2}_{*}\neq 0,\lambda^{2}_{*}\neq 0$, and two $O(3)$ QED-GNY fixed points with nonzero values for all couplings. The only stable fixed point among all eight is one of the latter two, with fixed-point couplings given by

[TABLE]

where we define

[TABLE]

These three coupling constants are positive for all values of $N_{f}$—including the case $N_{f}=1$, appropriate for spin-1/2 microscopic degrees of freedom—as in the chiral Ising QED3-GNY model Janssen and He (2017); Ihrig et al. (2018a); Zerf et al. (2018); thus the fixed point is physical and describes a valid QCP. This establishes that at one-loop order, the Néel-ASL transition is continuous and its critical properties can be systematically calculated in the $\epsilon$ expansion. From the four-loop order results in Eq. (58), a power series expansion of the critical exponents up to fourth order in $\epsilon$ can be obtained. The first exponent one can calculate is the stability exponent $\omega$, which controls leading corrections to scaling and would appear, for instance, in the subleading dependence of the AF susceptibility on the critical tuning parameter near the QCP. Mathematically, it corresponds to the smallest eigenvalue of the stability matrix, i.e., the Jacobian matrix $\mathcal{J}_{\alpha\alpha^{\prime}}=\partial\beta_{\alpha^{2}}/\partial\alpha^{\prime 2}$, $\alpha,\alpha^{\prime}=e,g,\lambda$ of derivatives of the beta functions evaluated at the fixed point (68)-(70). We find $\omega=\epsilon+\mathcal{O}(\epsilon^{2})$ at one-loop order. The full four-loop result for general $N_{f}$ is given in Ref. Sup ; here we present it only for $N_{f}=1$, appropriate to the Néel-ASL transition:

[TABLE]

In this and later $\mathcal{O}(\epsilon^{4})$ expressions for the $N_{f}=1$ critical exponents, coefficients are given to four significant digits.

III.4 Order parameter anomalous dimension

Microscopically, the order parameter anomalous dimension $\eta_{\phi}$ controls the long-distance behavior of AF correlations at the QCP,

[TABLE]

where $\bm{N}(r)\sim(-1)^{x+y}{\bm{S}}_{r}$ and $|r-r^{\prime}|$ is much greater than the lattice constant (but smaller than the system size, for a finite-size system approaching the thermodynamic limit). In a physical spin system this exponent would be extracted from an analysis of the width of the $(\pi,\pi)$ spin structure factor, measured in neutron scattering. In the field theory it is computed by evaluating $\gamma_{\phi}$ in Eq. (57) at the $O(3)$ QED-GNY fixed point (68)-(70),

[TABLE]

As for the beta functions, we express the anomalous dimensions as sums of contributions at fixed loop order,

[TABLE]

We find

[TABLE]

with the four-loop contribution $\gamma_{\phi}^{\text{(4L)}}$ given electronically Sup . At one-loop order, we find

[TABLE]

Evaluating the full four-loop result Sup for $N_{f}=1$, we obtain

[TABLE]

III.5 Correlation length exponent

In addition to the anomalous dimension $\gamma_{\phi}$ of the order parameter field $\bm{\phi}$, we also compute $\gamma_{m^{2}}\equiv\gamma_{\phi^{2}}-\gamma_{\phi}$ which appears in the beta function for the scalar mass squared,

[TABLE]

Evaluated at the QCP, the anomalous dimension $\gamma_{m^{2}}$ is related to the correlation length exponent $\nu$,

[TABLE]

which controls the divergence of the zero-temperature correlation length upon approach to the QCP. The contributions to $\gamma_{m^{2}}$ up to three-loop order are given by

[TABLE]

and the four-loop contribution $\gamma_{m^{2}}^{\text{(4L)}}$ is given in electronic format Sup . At one-loop order, we find

[TABLE]

while the full four-loop result Sup for $N_{f}=1$ is

[TABLE]

III.6 CDW and VBS susceptibility exponents

Besides the Néel order parameter field $\bm{\phi}$, which can be identified with the bilinear $\overline{\Psi}\bm{\sigma}\Psi$, other fermion bilinears will develop universal power-law correlations at the QCP. In the original $d=3$ theory with two-component spinors $\psi_{\alpha\sigma}$, these are Lorentz-invariant bilinears of the form $\overline{\psi}\tau_{i}\sigma_{j}\psi$ where the $\tau_{i}$ are Pauli matrices acting on the nodal/valley index $\alpha=\pm$. Besides the spin-triplet Néel bilinear $\overline{\Psi}\bm{\sigma}\Psi=\overline{\psi}\tau_{3}\bm{\sigma}\psi$ already considered, in this section we consider the three spin-singlet bilinears of the form $\overline{\psi}\tau_{i}\psi$, $i=1,2,3$. The $\overline{\psi}\psi$ spin-singlet bilinear will be discussed in Sec. IV.3.

The physical meaning of these bilinears can be ascertained from an analysis of their transformation properties under space-group and time-reversal symmetries via the projective symmetry group (PSG) Wen (2002) (see Appendix A). First, the transformation properties of the spin-triplet bilinear $\overline{\psi}\tau_{3}\bm{\sigma}\psi$ agree with those of the microscopic Néel order parameter and confirm the identification obtained in the naive continuum limit (Sec. II.1). The spin-singlet bilinear $\overline{\psi}\tau_{3}\psi=\overline{\Psi}\Psi$ transforms in the same way as the Néel bilinear under space-group symmetries but oppositely under time reversal; it thus corresponds to a staggered density or CDW order parameter for the lattice $c$ fermions,

[TABLE]

As discussed in Ref. Hermele et al. (2005), in the language of an underlying spin model this operator is symmetry equivalent to the staggered component of the density of skyrmions in the Néel field. The spin-singlet bilinears $\overline{\psi}\tau_{1}\psi$ and $\overline{\psi}\tau_{2}\psi$ transform into each other under $C_{4}$ rotations, are both even under time reversal, and transform oppositely under reflections and unit translations. In four-component notation with the choice of $4\times 4$ gamma matrices introduced in Sec. II.1, these can be written as $\overline{\psi}\tau_{1}\psi=i\overline{\Psi}\Gamma_{3}\Psi$ and $\overline{\psi}\tau_{2}\psi=-i\overline{\Psi}\Gamma_{5}\Psi$. These two bilinears can thus be identified with the two components of the VBS order parameter $\bm{V}\equiv(V_{x},V_{y})$, where $V_{x}\sim i\overline{\Psi}\Gamma_{5}\Psi$ and $V_{y}\sim i\overline{\Psi}\Gamma_{3}\Psi$ correspond to columnar VBS order in the $x$ and $y$ directions, respectively. In the microscopic theory (II.1), these can be defined as Xu et al. (2019)

[TABLE]

The microscopic CDW and VBS operators are functions of on-site bilinears of the form $c_{r\sigma}^{\dagger}c_{r\sigma^{\prime}}^{\phantom{\dagger}}$ and are thus manifestly gauge invariant.

From the operator identifications above we therefore expect that at the QCP, the CDW and VBS order parameters should develop universal power-law correlations at long distances $|r-r^{\prime}|\gg 1$ (measured in units of the lattice constant),

[TABLE]

where

[TABLE]

In the language of an underlying spin model, the VBS and CDW correlation functions would correspond to four-spin and six-spin correlation functions, respectively, which can in principle be measured by two-magnon and three-magnon inelastic light scattering, respectively Fleury (1969). Note that these correlations are already power-law in the ASL phase, as observed numerically for the Néel and VBS order parameters in Ref. Xu et al. (2019). However, the exponents in Eq. (92)-(93), as well as the exponent $\Delta_{\text{AF}}=\Delta_{\phi}$ for AF spin-spin correlations (i.e., the scaling dimension of the order parameter field $\bm{\phi}$), reflect the particular conformal field theory associated with the Néel-ASL QCP and are thus predicted to be different than those in the ASL phase, which corresponds to the pure conformal QED3 fixed point.

We further observe that while $i\overline{\Psi}\Gamma_{5}\Psi$ is a Lorentz (pseudo)scalar in $d=4$ dimensions, $i\overline{\Psi}\Gamma_{3}\Psi$ is not, but rather corresponds to one component of a Lorentz vector $i\overline{\Psi}\Gamma_{\mu}\Psi$. However, both are Lorentz scalars in $d=3$ dimensions. Therefore, the $\epsilon$ expansion will predict different scaling dimensions for $i\overline{\Psi}\Gamma_{5}\Psi$ and $i\overline{\Psi}\Gamma_{3}\Psi$, although they are expected to have the same scaling dimension in the physical $d=3$ theory since they transform into each other under the action of microscopic symmetries. For the purposes of determining $\Delta_{\text{VBS}}$, in the $\epsilon$ expansion we will therefore only calculate the dimension of the Lorentz pseudoscalar $i\overline{\Psi}\Gamma_{5}\Psi$, which does not have the unphysical property of transforming as a Lorentz vector in $d=4$ dimensions. We remind the reader that, as explained in the first paragraph of Sec. III, the entirety of Sec. III (and in fact the entire paper, except for Sec. II.2) is concerned with the Néel-ASL transition on the square lattice, for which the $\epsilon$-expansion is a meaningful approach. All the exponents are well defined in physical $d=2+1$ dimensions. We are adopting a conservative approach, where in the $\epsilon$-expansion analysis we only calculate the scaling dimension of quantities that remain Lorentz invariant when continuing from three to four spacetime dimensions. By contrast, in the $1/N_{f}$ expansion in $d=3$ the scaling dimensions of $i\overline{\Psi}\Gamma_{5}\Psi$ and $i\overline{\Psi}\Gamma_{3}\Psi$ should agree order by order in $1/N_{f}$. On the other hand, $\Delta_{\text{CDW}}$ can be meaningfully identified with $\Delta_{\overline{\Psi}\Psi}$ computed in both the $\epsilon$ and $1/N_{f}$ expansions with four-component spinors, since $\overline{\Psi}\Psi$ is a Lorentz scalar in both $d=4$ and $d=3$ dimensions.

The presence of $\Gamma_{5}$ can create issues in the $\epsilon$-expansion if one insists on preserving four-dimensional identities such as $\Gamma_{5}=\frac{1}{4!}\epsilon_{\mu\nu\lambda\sigma}\Gamma_{\mu}\Gamma_{\nu}\Gamma_{\lambda}\Gamma_{\sigma}$, which cannot be easily generalized to noninteger $d$ ’t Hooft and Veltman (1972). In our case, we are interested in describing a quantum critical point in $d=3$, for which the intrinsically four-dimensional tensor $\epsilon_{\mu\nu\lambda\sigma}$ does not exist. It is therefore sufficient to use the naive prescription according to which $\{\Gamma_{5},\Gamma_{\mu}\}=0$ in all $d$ dimensions Chanowitz et al. (1979). In this prescription evanescent operators Dugan and Grinstein (1991) are not generated when renormalizing the $i\overline{\Psi}\Gamma_{5}\Psi$ operator. Evanescent operators can be generated when renormalizing four-fermion operators Di Pietro and Stamou (2018), but we will not concern ourselves with the renormalization of such operators here.

We now turn to the computation of the scaling dimensions $\Delta_{\overline{\Psi}\mathcal{M}\Psi}$ of fermion bilinears $\overline{\Psi}\mathcal{M}\Psi$ (here $\mathcal{M}$ is the identity matrix or $i\Gamma_{5}$). As discussed above, such scaling dimensions describe the universal power-law decay in real space of various equal-time correlation functions (i.e., thermodynamic susceptibilities) at the quantum critical point. The way those scaling dimensions are calculated here is similar to how a linear-response susceptibility is usually defined/measured in statistical mechanics, i.e., one turns on a weak perturbation and computes the response to this perturbation, in the limit that the strength of the perturbation goes to zero. Here, to compute the fermion bilinear scaling dimensions, one must add a small, perturbative fermion mass term (really, turn on a weak staggered or VBS order), and compute the ratio of the induced order to the applied perturbation. (Since the added mass is infinitesimally weak, the fermions remain Dirac fermions, albeit slightly massive.) As usual with linear response, however, the ratio itself is a property of the system in the absence of the perturbation. Here, the susceptibility exponents are universal properties of the quantum critical point itself, where there are no mass terms and the fermions are gapless. To put these ideas into practice, we follow the method discussed in Ref. Zerf et al. (2018), whereby the bilinear is added to the bare (renormalized) Lagrangian with a coefficient $M_{0}$ ($M$), the two coefficients being related by $M=M_{0}\mu^{-1}Z_{M}^{-1}Z_{\Psi}$. The renormalized mass $M$ obeys the RG equation

[TABLE]

where $\gamma_{M}=\mu(d\ln Z_{M}/d\mu)$. Since $M$ is the coefficient of a gauge-invariant operator, the combination $\gamma_{M}-\gamma_{\Psi}$ must be gauge invariant, even though $\gamma_{M}$ and $\gamma_{\Psi}$ individually are not. The scaling dimension of the bilinear is then given by

[TABLE]

where $\gamma_{\overline{\Psi}\mathcal{M}\Psi}\equiv\gamma_{M}-\gamma_{\Psi}$. By contrast with the chiral Ising QED-GNY model Zerf et al. (2018), in which a $\overline{\Psi}\Psi$ insertion (i.e., adding this term to the Lagrangian) breaks a $\mathbb{Z}_{2}$ chiral symmetry and radiatively induces a relevant $\phi^{3}$ interaction in the renormalized Lagrangian, here the continuous $O(3)$ symmetry ensures that no additional operators which might mix with $\overline{\Psi}\mathcal{M}\Psi$ are radiatively induced by an insertion of this operator 222A radiatively induced interaction is one which is not present in the original Lagrangian but appears at higher orders in perturbation theory due to loop diagrams..

The CDW susceptibility exponent is given in terms of the normal mass operator anomalous dimension $\gamma_{\overline{\Psi}\Psi}$, given up to three-loop order by

[TABLE]

The four-loop term $\gamma_{\overline{\Psi}\Psi}^{\text{(4L)}}$ is given electronically Sup . Contributions are found to be independent of the gauge-fixing parameter $\xi$ at each loop order, which constitutes a strong check on the calculation. Evaluating Eq. (95) at the QCP (68)-(70), we obtain at one-loop order

[TABLE]

with the full four-loop expression given in Ref. Sup . For $N_{f}=1$, we obtain

[TABLE]

Similarly, the VBS susceptibility exponent is obtained from the axial mass operator anomalous dimension $\gamma_{i\overline{\Psi}\Gamma_{5}\Psi}$. Up to three-loop order we obtain

[TABLE]

with the four-loop contribution $\gamma_{i\overline{\Psi}\Gamma_{5}\Psi}^{\text{(4L)}}$ given in electronic format Sup . At one-loop order for general $N_{f}$, we find

[TABLE]

while for $N_{f}=1$ at four-loop order, we obtain

[TABLE]

The full four-loop expression for general $N_{f}$ is given in Ref. Sup .

Finally, for $N_{f}>1$, in addition to the flavor-singlet bilinears $\overline{\Psi}\Psi$ and $i\overline{\Psi}\Gamma_{5}\Psi$ one can define an $SU(N_{f})$ flavor-adjoint bilinear,

[TABLE]

where $T_{A}$ are the Hermitian generators of $SU(N_{f})$. This operator does not exist for $N_{f}=1$ and thus has no physical meaning in the context of the Néel-ASL transition, but the computation of its scaling dimension for $N_{f}>1$ is an interesting problem in quantum field theory. The corresponding anomalous dimension $\gamma_{\overline{\Psi}T_{A}\Psi}$ is given up to four-loop order in Ref. Sup . At one-loop order, $\Delta_{\overline{\Psi}T_{A}\Psi}$ is equal to the flavor-singlet, normal mass bilinear dimension (99). A difference in scaling dimensions between singlet and adjoint bilinears begins appearing only at four-loop order,

[TABLE]

As discussed in Ref. Zerf et al. (2018), the difference between singlet and adjoint mass dimensions comes from bilinear insertions in closed fermion loops appearing in the loop expansion of the fermion two-point function; such closed loops vanish for adjoint insertions since the $SU(N_{f})$ generators are traceless, but are generally nonzero for singlet insertions. However, here the tracelessness of the spin $SU(2)$ generators appearing at each Yukawa vertex implies that nonvanishing singlet insertions in closed fermion loops occur at higher loop order than in the chiral Ising QED-GNY model, since such closed loops must contain an even number of Yukawa vertices to not vanish in the chiral Heisenberg case.

III.7 Fermion bilinear dimensions in the chiral Heisenberg GNY model

Turning off the gauge coupling $e^{2}$ in the Lagrangian (III.1), we obtain the pure chiral Heisenberg GNY model, which for $N_{f}=1$ describes the universal critical properties of the semimetal-AF insulator quantum phase transition in the half-filled Hubbard model on the honeycomb and $\pi$-flux square lattices Sorella and Tosatti (1992); Herbut (2006); Honerkamp (2008); Herbut et al. (2009); Sorella et al. (2012); Ulybyshev et al. (2013); Assaad and Herbut (2013); Janssen and Herbut (2014); Parisen Toldin et al. (2015); Otsuka et al. (2016); Buividovich et al. (2018). In this limit, the Dirac fermions can be interpreted as real electronic excitations in a semimetal, rather than fractionalized spinon excitations in a spin system. Whereas the usual exponents ($\nu$, $\eta_{\phi}$, and the fermion anomalous dimension $\eta_{\Psi}$) for the chiral Heisenberg GNY model have been computed both in the $\epsilon$- Rosenstein et al. (1993); Zerf et al. (2017) and large-$N_{f}$ Gracey (2018c) expansions, to our knowledge scaling dimensions of fermion bilinears have thus far not been computed. In particular, the $\overline{\Psi}\Psi$ bilinear corresponds in the long-wavelength limit to the CDW order parameter, i.e., Eq. (88) on the $\pi$-flux square lattice, and the Semenoff mass Semenoff (1984) on the honeycomb lattice. Correlations of this bilinear can in principle be computed in sign-problem-free QMC simulations of the half-filled repulsive Hubbard model on the aforementioned lattices, and at the semimetal-AF QCP should decay at long distances according to Eq. (90) with an exponent $\Delta_{\text{CDW}}=\Delta_{\overline{\Psi}\Psi}^{\text{cH-GNY}}$ given by

[TABLE]

For the chiral Heisenberg GNY model with general $N_{f}$, we obtain at one-loop order

[TABLE]

A four-loop expression for general $N_{f}$ is given electronically in Ref. Sup .

Similarly, the Hermitian $i\overline{\Psi}\Gamma_{3}\Psi$ and $i\overline{\Psi}\Gamma_{5}\Psi$ bilinears can be viewed as the two components of the columnar VBS order parameter on the $\pi$-flux square lattice, or alternatively as the real and imaginary parts of a $\mathbb{Z}_{3}$-symmetric complex order parameter describing Kekulé VBS order on the honeycomb lattice Hou et al. (2007); Ryu et al. (2009); Roy and Herbut (2010). At the semimetal-AF QCP of the half-filled repulsive Hubbard model on those lattices the correlation function of these operators should decay as in Eq. (91), with the exponent $\Delta_{\text{VBS}}=\Delta_{i\overline{\Psi}\Gamma_{5}\Psi}^{\text{cH-GNY}}$ given by

[TABLE]

For general $N_{f}$, the one-loop result is

[TABLE]

and the full four-loop result is given in electronic format Sup . By contrast with the corresponding exponents at the Néel-ASL QCP, which would describe multi-spin correlations in a microscopic spin system, at the semimetal-AF insulator QCP these exponents could be accessed by measuring ordinary two-point correlations of one-body operators for electrons, e.g., density-density and current-current correlations.

Turning finally to the $SU(N_{f})$ flavor-adjoint bilinear (106), the four-loop diagrams responsible for the difference (95) in scaling dimensions between the singlet and adjoint bilinears are $\mathcal{O}(e^{2}g^{6})$ and vanish for the pure chiral Heisenberg GNY model with $e^{2}=0$. Thus for the latter model this difference in scaling dimensions can at most be $\mathcal{O}(\epsilon^{5})$.

IV Large-$N_{f}$ expansion

To complement the $\epsilon$-expansion studies just discussed, we now turn to the large-$N_{f}$ expansion (for a recent review, see Ref. Gracey (2018d)). In contrast with Sec. III, where $\epsilon=4-d$ is treated as a formal expansion parameter but $N_{f}$ is arbitrary, here we work in fixed $d$ dimensions but treat $1/N_{f}$ as a formal expansion parameter. That is, we consider a system with $N_{f}$ flavors of $SU(2)$ doublets of four-component Dirac fermions (i.e., a total of $2N_{f}$ four-component Dirac fermions), where we assume $N_{f}\gg 1$, and perform perturbation theory in powers of $1/N_{f}$; at the end of the calculation we then consider the physically appropriate limit $N_{f}\rightarrow 1$. (Larger values of $N_{f}$ may be physically accessible in systems of large-spin ultracold fermionic atoms such as 87Sr Gorshkov et al. (2010); DeSalvo et al. (2010).) The large-$N_{\!f}$ formalism we use has been detailed in Refs. Gracey (1992, 1993b); *gracey2018b; Gracey (1993c, 2018b) for the chiral Ising QED-GN(Y) model. We refer the reader to those articles for more details as well as to the original papers of Vasil’ev et al. Vasil’ev et al. (1981a, b), where the large-$N_{\!f}$ critical point approach was developed for scalar theories with $N_{f}$ flavors in $d$ dimensions. In essence, the zeroth-order ($N_{f}=\infty$) approximation corresponds to calculating the gauge-field and scalar-field propagators in the random-phase approximation (RPA), and corrections at higher orders in $1/N_{f}$ come from diagrams containing self-energy and vertex corrections computed with the RPA-resummed scalar-field and gauge-field propagators. Thus, the diagrams are grouped according to their powers in $1/N_{f}$, irrespective of their loop order, which is the classification used for grouping diagrams in the $\epsilon$ expansion. The elegance of the large-$N_{f}$ method is such that one carries out all computations directly at the $d$-dimensional fixed point where there is scaling behavior. The core Lagrangian for determining large-$N_{\!f}$ critical exponents is that which governs the universality class of the $O(3)$ QED-GNY fixed point of the four-dimensional Lagrangian (III.1) used for the perturbative computations. The essence of that Lagrangian is that it involves the kinetic terms of the fields as well as the $3$-point interactions between the matter and force fields. By the latter we include not only the photon but also the scalar field $\phi^{a}$. For our case there are only two such interactions and the Euclidean Lagrangian of the universal theory used for the large-$N_{\!f}$ analysis is

[TABLE]

where sums over repeated indices $i=1,\ldots,N_{f}$ and $a=1,2,3$ are understood. We have not included any additional terms since there would be an infinite number of operators built from the three fields. Only a few of these operators would be relevant in a fixed dimension; in particular, the kinetic terms for the scalar and gauge fields in Eq. (III.1) are irrelevant at criticality in the large-$N_{f}$ limit. Also the coupling constants have been rescaled out of the interactions to produce the new fields $\tilde{\phi}^{a}$ and $\tilde{A}_{\mu}$ since the formalism applies at criticality and there the couplings do not run Vasil’ev et al. (1981a, b). Excluding one or the other of the two interactions in Eq. (112) would lead to critical exponents of a different universality class. At criticality the canonical dimensions of the fields are fixed by the dimensionlessness of the action in $d$ dimensions and this leads to the asymptotic scaling forms of the propagators Gracey (1992, 2018b),

[TABLE]

in the approach to the fixed point in coordinate space. Here $A$, $B_{A}$ and $B_{\phi}$ are coordinate-independent amplitudes and the full dimension of each field is defined by

[TABLE]

where $\chi_{A}$ and $\chi_{\phi}$ are the anomalous dimensions of the respective vertices of (112) and we retain the definition $d$ $=$ $2\mu$ of Refs. Gracey (2018b); Vasil’ev et al. (1981a) for consistency with other large-$N_{\!f}$ work using this formalism. To avoid confusion with the notation for the $\beta$-functions we use a hat for the full dimension of the field exponents. The vertex dimensions as well as the other exponents such as $\eta$ are functions of $d$ and $N_{\!f}$ and each exponent can be expanded in a series of the form

[TABLE]

Moreover these critical exponents are related to the $\epsilon$-expansion of their respective RG functions at the $O(3)$ QED-GNY fixed point found in Sec. III.3. As such this provides a non-trivial check on the four-loop renormalization. Since the large-$N_{\!f}$ formalism of Refs. Vasil’ev et al. (1981a, b) operates purely at the critical point all our $1/N_{\!f}$ results will be exclusively in the Landau gauge. The gauge parameter, $\xi$, of our linear gauge fixing can be regarded as a second coupling constant and at criticality its associated RG function, which is in effect a $\beta$-function, defines a critical coupling. In this case it corresponds to the Landau gauge.

IV.1 Critical exponents

Using this formalism we have extended the $\mathcal{O}(1/N_{\!f}^{2})$ computations of Ref. Gracey (2018b) to the case where the Pauli matrices appear in the scalar-fermion vertex. Since Refs. Gracey (1993c, 2018b) already record the values for all the large-$N_{\!f}$ master integrals contributing to the scaling forms of the various Green’s functions needed to find the operator dimension at $\mathcal{O}(1/N_{\!f}^{2})$, we merely record the outcome of this straightforward exercise. The key building block is the fermion anomalous dimension, which at leading order is

[TABLE]

in Landau gauge, where $\Gamma(z)$ is the gamma function. Having established this result, the two leading-order vertex dimensions are required for the solution of the $\mathcal{O}(1/N_{\!f}^{2})$ Schwinger-Dyson equations. We find

[TABLE]

The simplicity of the latter equation is a reflection that in the critical large-$N_{\!f}$ setup the photon has a full dimension of unity which correctly shows the consistency of the Ward-Takahashi identity. In the RG approach this is equivalent to the all-orders statement that $\gamma_{A}=\epsilon$ at any fixed point with $e_{*}^{2}\neq 0$, which follows from Eq. (53). In addition we have computed $1/\nu$ to $\mathcal{O}(1/N_{\!f})$ and found

[TABLE]

With the leading order exponents for $\eta$, $\chi_{A}$, and $\chi_{\phi}$ we have evaluated several exponents to the next order finding

[TABLE]

where $\psi(z)=\Gamma^{\prime}(z)/\Gamma(z)$ is the Euler digamma function and we define

[TABLE]

As these expressions are valid in $d$ dimensions, expanding them in powers of $\epsilon$ with $d=4-\epsilon$ and comparing with the corresponding RG functions at criticality we find full agreement, which is a solid check on both the perturbative and large-$N_{\!f}$ computations. The fermion anomalous dimension $\eta$ obtained from Eqs. (116) and (119) agrees with the Landau-gauge ($\xi=0$) fermion anomalous dimension obtained from the four-loop RG function $\gamma_{\Psi}$ defined in Eq. (57). The exponent $\hat{\beta}_{\phi}$ corresponds to the full scaling dimension of the scalar field, which is given by $(d-2+\eta_{\phi})/2$ in $d$ dimensions where the anomalous dimension $\eta_{\phi}$ corresponds to the RG function $\gamma_{\phi}$ at criticality [Eq. (74)]. Finally, Eq. (118) agrees with the inverse correlation length exponent found in Eq. (82). In $d=3$ dimensions we obtain

[TABLE]

or, numerically

[TABLE]

using $\eta_{\phi}=2\left[2-\mu-\left(\eta+\chi_{\phi}\right)\right]$.

IV.2 CDW and VBS susceptibility exponents

The final exercise is to determine the scaling dimension of fermion bilinears in the large-$N_{f}$ expansion. We begin with the flavor-singlet bilinear $\overline{\Psi}\Psi$, which determines the CDW susceptibility exponent [see Eqs. (90) and (92)], as well as the adjoint bilinear $\overline{\Psi}T_{A}\Psi$. We shall compute these at $\mathcal{O}(1/N_{\!f}^{2})$, which requires the expression for $\eta_{2}$. In this case there is a subtle difference between the setup for Eq. (112) and that of the theory where the Pauli matrix of the scalar-fermion vertex is absent. The latter was considered in Ref. Gracey (2018b) and in that case the $\mathcal{O}(1/N_{\!f}^{2})$ expressions for the adjoint and singlet operator dimensions were different. This was because of the graphs of Figure $7$ of Ref. Gracey (2018b). These correspond to the insertion of the operator in a closed fermion loop. For the adjoint case these graphs are absent due to the trace over the matrix of the operator. The operator of the singlet case in Ref. Gracey (2018b) is the fermion mass operator and the corresponding graphs contribute leading to different $\mathcal{O}(1/N_{\!f}^{2})$ exponents. However, in the case of Eq. (112) the presence of the Pauli matrix in the scalar-fermion vertex means that there is no contribution from the corresponding graphs of Figure $7$ of Ref. Gracey (2018b). This does not mean that the singlet and adjoint operator dimensions are equivalent in the perturbative case. In fact they are different, with the discrepancy apparent at $\mathcal{O}(1/N_{\!f}^{3})$, as can be proved explicitly by expanding Eq. (107) in inverse powers of $N_{f}$. Computations at $\mathcal{O}(1/N_{\!f}^{3})$ are however beyond the reach of the current large-$N_{\!f}$ formalism. The outcome is that we have the leading-order exponent

[TABLE]

and

[TABLE]

at next-to-leading order, for both singlet and adjoint bilinears. We have checked that the $\epsilon$-expansion of both are in full agreement with the four-loop order results obtained from Eq. (95). The next stage is to deduce the three-dimensional values

[TABLE]

or

[TABLE]

where

[TABLE]

In contrast to the corresponding expression in Ref. Gracey (2018b) the coefficient of the $\mathcal{O}(1/N_{\!f}^{2})$ term is significantly smaller than that at $\mathcal{O}(1/N_{\!f})$. This may mean that exponent estimates for relatively low values of $N_{\!f}$ could be reasonably reliable. However, this should not be interpreted as a general statement regarding the relative reliability of the large-$N_{f}$ expansion at small $N_{f}$ for the Heisenberg versus Ising theories, since the coefficients of $\mathcal{O}(1/N_{f}^{3})$ and higher-order terms remain unknown in these theories and thus cannot be compared.

We next turn to the axial mass operator $i\overline{\Psi}\Gamma_{5}\Psi$, which controls the decay of VBS correlations [see Eqs. (91) and (93)] at the QCP. Expressing its scaling dimension as $\Delta_{i\overline{\Psi}\Gamma_{5}\Psi}=2\mu-1+\eta+\eta_{\mathcal{O}_{5}}$, the leading $\mathcal{O}(1/N_{f})$ correction is

[TABLE]

with the next-to-leading order term given by

[TABLE]

In $d=3$ dimensions, this evaluates to

[TABLE]

If the result for arbitrary $d$ is expanded in powers of $\epsilon$, where $d=4-\epsilon$, then the series that ensues agrees with the large-$N_{f}$ expansion of the four-loop order result in Eq. (95).

IV.3 QAH susceptibility exponent and Aslamazov-Larkin diagrams

We now turn to the spin-singlet bilinear $\overline{\psi}\psi=i\overline{\Psi}\Gamma_{3}\Gamma_{5}\Psi$, which is odd under time reversal $\mathcal{T}$ and reflection $R_{x}$ (i.e., parity in 2+1 dimensions) as shown in Appendix A. Because of the presence of the $\Gamma_{3}$ matrix, this bilinear is not Lorentz invariant in $d=4-\epsilon$ dimensions but is Lorentz invariant in strict $d=3$ dimensions, where the spacetime index $\mu=0,1,2$ only. As a result, we do not calculate its dimension in the $\epsilon$ expansion, but can compute it in the large-$N_{f}$ expansion in fixed $d=3$, using the relation $\Gamma_{5}=\Gamma_{0}\Gamma_{1}\Gamma_{2}\Gamma_{3}$. Physically, this operator corresponds to the continuum limit of a QAH or Haldane mass term Haldane (1988), first discussed in the context of the $\pi$-flux phase on the square lattice by Wen, Wilczek, and Zee Wen et al. (1989):

[TABLE]

The pure imaginary, diagonal next-nearest-neighbor hopping terms break time-reversal symmetry and give a mass to the Dirac fermions. We define two Wilson line operators as

[TABLE]

which obey $W_{r,r+{\bm{a}}_{1}}^{\dagger}=W_{r+{\bm{a}}_{1},r}$ and $W_{r,r+{\bm{a}}_{2}}^{\dagger}=W_{r+{\bm{a}}_{2},r}$. Such Wilson line insertions Hermele et al. (2005) are required to ensure that $\mathcal{O}_{\text{QAH}}(r)$ is invariant under lattice $U(1)$ gauge transformations $c_{r\sigma}\rightarrow e^{i\varphi_{r}}c_{r\sigma}$, $\theta_{r,r+\hat{\mu}}\rightarrow\theta_{r,r+\hat{\mu}}+\varphi_{r}-\varphi_{r+\hat{\mu}}$, and their particular form has been chosen to preserve all lattice symmetries: one can check that Eq. (IV.3) transforms precisely as the continuum $\overline{\psi}\psi$ bilinear under the microscopic PSG transformations in Eq. (A.11). As result we expect that at the Néel-ASL QCP, correlations of the QAH operator (IV.3) as computed in QMC would decay in power-law fashion at long distances,

[TABLE]

with

[TABLE]

where $\Delta_{i\overline{\Psi}\Gamma_{3}\Gamma_{5}\Psi}$ is the dimension of the continuum bilinear.

In the preceding subsections critical exponents have been computed in arbitrary $d$-dimensional spacetime utilizing the fact that the trace of the product $\Gamma_{\mu_{1}}\cdots\Gamma_{\mu_{2n+1}}$ of an odd number of $4\times 4$ gamma matrices is zero. In fixed $d=3$ dimensions with two-component spinors, this assumption is invalid since $\mathop{\mathrm{tr}}\gamma_{\mu}\gamma_{\nu}\gamma_{\lambda}\propto\epsilon_{\mu\nu\lambda}$, which was shown to give additional contributions to the large-$N_{f}$ critical exponents of the chiral Ising QED3-GNY model from Aslamazov-Larkin diagrams Boyack et al. (2019); Benvenuti and Khachatryan (2019). In the four-component formulation with $d=3$, we expect Aslamazov-Larkin diagrams to again contribute in the presence of a $\Gamma_{3}\Gamma_{5}$ insertion, since $\mathop{\mathrm{tr}}\Gamma_{3}\Gamma_{5}\Gamma_{\mu}\Gamma_{\nu}\Gamma_{\lambda}\propto\epsilon_{\mu\nu\lambda}$. The calculation of the Aslamazov-Larkin diagrams for the chiral-Ising case was performed in Ref. Boyack et al. (2019), and so those results can be adapted to the present case by incorporating the additional $SU(2)$ spin structure present in the Feynman rules. As a result of the modified Feynman rules, the large-$N_{f}$ scalar propagator is twice the value of that appearing in Ref. Boyack et al. (2019), while the photon propagator is now one half of its value in the chiral-Ising case. In reference to the Feynman diagrams in Ref. Boyack et al. (2019), we note that since there are two spin components, the diagram in Fig. 4 is multiplied by two, whereas the spin factor and the modified photon propagator factor cancel one another to render Figs. 5(a-e) equal to their original value. For the diagrams in Figs. 6(a-c), the Pauli matrices in the Yukawa vertex now lead to an additional factor of $\frac{1}{4}\mathrm{tr}\sigma^{a}\sigma^{b}\delta^{ab}=3/2$, and thus the diagrams as a whole are $2\times 3/2$ larger than their original values. The diagrams in Figs. 7(a,d) still give zero total contribution. As a result, the QAH susceptibility exponent in Eq. (135) is given to $\mathcal{O}(1/N_{f})$ by

[TABLE]

IV.4 Fermion bilinear dimensions in the chiral Heisenberg GN(Y) model

While the fermion ($\eta$) and boson ($\eta_{\phi}$) anomalous dimensions as well as the inverse correlation length exponent $1/\nu$ have been computed in the large-$N_{f}$ expansion for the pure chiral Heisenberg GN(Y) model Gracey (2018c), the scaling dimensions of fermion bilinears have thus far not been computed using this method. Using the same formalism as in previous subsections, we have computed the dimension of the flavor-singlet $\overline{\Psi}\Psi$ and flavor-adjoint $\overline{\Psi}T_{A}\Psi$ bilinears at $\mathcal{O}(1/N_{f}^{2})$, which are the same at this order, as well as the dimension of the axial bilinear $i\overline{\Psi}\Gamma_{5}\Psi$ at $\mathcal{O}(1/N_{f}^{2})$. The general expressions are again Eq. (128) and its axial counterpart, together with Eq. (115), but where the exponents $\eta$, $\eta_{\mathcal{O}}$, and $\eta_{\mathcal{O}_{5}}$ must be computed at the critical fixed point of the pure chiral Heisenberg GN(Y) theory. The fermion exponents $\eta_{1}$ and $\eta_{2}$ were computed previously Gracey (2018c), but are reproduced here for convenience and accounting for the change in convention regarding the definition of flavor number ($N$ vs $N_{f}$). We obtain

[TABLE]

at $\mathcal{O}(1/N_{f})$, and

[TABLE]

at $\mathcal{O}(1/N_{f}^{2})$, defining

[TABLE]

In $d=4-\epsilon$ dimensions the results for $\Delta_{\overline{\Psi}\Psi}^{\text{cH-GNY}}$ and $\Delta_{i\overline{\Psi}\Gamma_{5}\Psi}^{\text{cH-GNY}}$ thus found agree with the corresponding four-loop results in Sec. III.7. The latter are obtained from Eq. (95), but evaluated at the pure GNY fixed point. For $d=3$, we find

[TABLE]

As discussed in Sec. III.7, for $N_{f}=1$ those exponents control the decay of CDW and VBS correlations at the semimetal-AF QCP in the repulsive Hubbard model on the honeycomb and $\pi$-flux square lattices, and are in principle accessible to sign-problem-free QMC simulations.

IV.5 Monopole operators

The microscopic theory studied in Ref. Xu et al. (2019) is defined on a lattice and thus has a compact $U(1)$ gauge group, as a result of the magnetic flux being defined only modulo $2\pi$ on the lattice. In such theories total magnetic flux is not conserved, unlike in noncompact $U(1)$ gauge theories, and monopole operators—gauge-invariant local operators which insert $2\pi q$ magnetic flux at a point in space, with $q\in\mathbb{Z}$—generically appear in the effective low-energy Lagrangian Polyakov (1975, 1977, 1987), provided they are invariant under all microscopic symmetries. If allowed monopole operators are relevant under the RG, they will proliferate, gap out the $U(1)$ gauge field, and confine excitations with nonzero gauge charge (here, the fermionic spinons); if irrelevant, compactness effects disappear at long distances. In a phase where the fermions are gapped by a symmetry-breaking order parameter, e.g., the Néel phase in our case, monopoles are known to be relevant and do proliferate Polyakov (1975, 1977, 1987), leading to a conventional AF phase with no propagating fractionalized excitations (spinons or low-energy gauge fields). When the fermions are gapless, such as in the ASL phase or at the Néel-ASL QCP, the monopoles are “dressed” by these gapless degrees of freedom and their scaling dimensions change. Those scaling dimensions can be computed in the large-$N_{f}$ expansion using the state-operator correspondence of conformal field theory Borokhov et al. (2002); a $2\pi q$ monopole operator is irrelevant if its scaling dimension $\Delta_{2\pi|q|}$ is greater than 3.

In the ASL phase, described by conformal QED3, in the large-$N_{f}$ limit the scaling dimension of the smallest ($q=\pm 1$) monopole operator allowed by the symmetries of the square lattice is $\Delta_{2\pi}^{\text{ASL}}=1.06N_{f}-0.0383+\mathcal{O}(1/N_{f})$ Pufu (2014), where $N_{f}$ is, according to our convention, the number of $SU(2)$ doublets of four-component Dirac spinors. (Monopoles with larger $|q|$ have larger scaling dimensions, at least at leading order in the $1/N_{f}$ expansion.) Considering only the first two terms in the large-$N_{f}$ expansion would suggest that monopoles are relevant for $N_{f}=1$, which would destabilize the ASL phase. However, the QMC results Xu et al. (2019); Wang et al. (2019) find a stable ASL phase with power-law correlations in accordance with the predictions of the noncompact theory. There are thus two possibilities: (1) the $\mathcal{O}(1/N_{f})$ and higher-order terms in $\Delta_{2\pi}^{\text{ASL}}$ become important at small $N_{f}$ and cause the monopoles to be irrelevant for $N_{f}=1$; (2) the monopoles are formally relevant, but the bare monopole fugacity $g_{0}$ is such that the length scale $L^{*}\sim ag_{0}^{-1/(3-\Delta_{2\pi}^{\text{ASL}})}$ beyond which monopole proliferation is observable (with $a$ the lattice constant) is much greater than the system sizes currently accessible in QMC simulations. (Equivalently, in the thermodynamic limit $L\rightarrow\infty$ monopole proliferation would occur only below a very small crossover temperature $T^{*}\sim c/L^{*}$ with $c$ some characteristic velocity.) In the latter scenario, critical scaling characteristic of the ASL phase would appear in an extended crossover regime $a\ll L\ll L^{*}$.

Turning to the Néel-ASL QCP itself, scaling dimensions of monopole operators at this critical point can again in principle be computed in the large-$N_{f}$ expansion. Such a calculation was recently performed Dupuis et al. (2019) for the theory in Eq. (34), describing the transition from the ASL to the ${\bm{q}}=0$ noncollinear state on the kagomé lattice. However, as already discussed, this theory is distinct from the theory (26) of the Néel-ASL transition on the square lattice. Generalizing to $N_{f}$ numbers of $SU(2)$ doublets of four-component Dirac fermions, which is required for a large-$N_{f}$ computation of monopole scaling dimensions, the theory (26) has an $SU(2)_{\text{spin}}\times SU(N_{f})_{\text{flavor}}\times U(1)_{\text{nodal}}$ symmetry, while the theory (34) has a larger $SU(2)_{\text{spin}}\times SU(2N_{f})_{\text{flavor-nodal}}$ symmetry. (On both lattices, and with our definition of $N_{f}$, the ASL phase has an enhanced $SU(4N_{f})$ symmetry.) That is because, as discussed in Sec. II.2, the bilinear that condenses on the kagomé lattice is proportional to the identity in nodal space, while the bilinear that condenses on the square lattice contains a $\tau_{3}$ nodal Pauli matrix and is only invariant under $U(1)$ nodal rotations about that direction. As a result, the large-$N_{f}$ monopole scaling dimensions computed for theory (34) in Ref. Dupuis et al. (2019) cannot be directly transposed to the Néel-ASL transition, since the two critical theories have different flavor symmetries. On the other hand, the QMC results Xu et al. (2019) suggest that a continuous Néel-ASL transition occurs at $N_{f}=1$. Indeed, Ref. Xu et al. (2019) calculated the flux energy per plaquette, for various values of the coupling strength $J$. This was observed to be a continuous quantity, which moreover was reproducible with increasing system size, and the absence of any singularity at the critical $J_{c}$ suggests a continuous phase transition. While one expects from the state-operator correspondence that $\Delta_{2\pi}^{\text{N\'{e}el-ASL}}\propto N_{f}$ in the large-$N_{f}$ limit, plus subleading corrections in $1/N_{f}$, two scenarios consistent with the QMC results can be envisaged as above, with the added uncertainty resulting from the unknown proportionality coefficient. (1) $\Delta_{2\pi}^{\text{N\'{e}el-ASL}}<3$ at $N_{f}=1$, and the transition becomes first order in the thermodynamic limit (assuming the ASL phase is itself stable) but the length scale $L^{**}\sim ag_{0}^{-1/(3-\Delta_{2\pi}^{\text{N\'{e}el-ASL}})}$ beyond which this happens is much greater than the system size $L$; critical behavior with the exponents computed here is observed for $a\ll L\ll L^{**}$. (2) $\Delta_{2\pi}^{\text{N\'{e}el-ASL}}>3$ at $N_{f}=1$ and true critical behavior persists in the thermodynamic limit.

In summary, the ultimate fate of monopoles in the deep infrared limit for $N_{f}=1$ remains unknown, both in the ASL phase and at the Néel-ASL QCP. The critical exponents computed in this work describe either true critical behavior in the thermodynamic limit, or quasicritical behavior below crossover lengths $L^{*}$ and $L^{**}$ which are larger than system sizes currently accessible in numerics.

V Resummed critical exponents

Mathematically, the $\epsilon$-expansion and the large-$N_{f}$ expansion both produce asymptotic series, with a radius of convergence of zero. Resummation procedures, which are consistent ways of assigning finite values to formally divergent series, are thus necessary to extract physical results from the finite-order expressions derived analytically in Sections III and IV. We now apply two particular approximate resummation methods to obtain estimates of critical exponents in fixed $d=3$ for general $N_{f}\geq 1$. A standard method to extrapolate finite-order $\epsilon$-expansion results to the physical dimension $d=3$ is the use of Padé approximants Kleinert and Schulte-Frohlinde (2001). For a given loop order $L$, the (one-sided) Padé approximants are defined by

[TABLE]

where $m$ and $n$ are two positive integers obeying $m+n=L$. The coefficients $a_{i}$ and $b_{j}$ are determined such that expanding the above function in powers of $\epsilon$ to $\mathcal{O}(\epsilon^{L})$ reproduces the $\epsilon$-expansion results. We also use the Padé-Borel method Kleinert and Schulte-Frohlinde (2001), which provides an alternative method to resummation of the $\epsilon$-expansion results. For an expansion given by $\Delta(\epsilon)=\sum_{k=0}^{\infty}\Delta_{k}\epsilon^{k}$, the Borel sum is defined as $B_{\Delta}(\epsilon)=\sum_{k=0}^{\infty}\frac{1}{k!}\Delta_{k}\epsilon^{k}$. The Borel transform is the following exact result:

[TABLE]

Since the $\epsilon$-expansion coefficients are known only to fourth order, in the above expression $B_{\Delta}$ is replaced by a Padé approximant at a given order. In the following we present the results of Padé extrapolation of the $\epsilon$-expansion for general $N_{f}$, as well as the corresponding Padé-Borel results for $N_{f}=1$ only. We also apply those two methods to the resummation of the $1/N_{f}$-expansion series, as done for instance in Ref. Gracey (2018b), for the case of $N_{f}=1$ only. In regard to the appropriateness of applying resummation techniques, we note that the large-$N_{f}$ expansion already at $\mathcal{O}(1/N_{f})$ agrees quite well with numerical QMC results for the ASL phase with $N_{f}=1$ in our notation (see Fig. 7 in Ref. Xu et al. (2019), which is $N_{f}=2$ in their notation). Since the ASL phase with $N_{f}=1$ is already well captured by the large-$N_{f}$ expansion, it is physically reasonable to investigate applying this method to phase transitions out of the ASL phase.

V.1 Chiral Heisenberg QED3-GNY model

The Padé approximants for $\eta_{\phi}$, $1/\nu$, and $\omega$ at two-loop (blue), three-loop (green), and four-loop (red) orders are shown in Figs. 2, 3, and 4 respectively. Only those Padé approximants that are nonsingular in the extrapolation region $0<\epsilon<1$ are shown in the figures. With increasing $N_{f}$ the various approximants generally converge towards each other as well as towards the $d=3$ large-$N_{f}$ result (black line). Furthermore, in Fig. 3 the spread of values predicted by the approximants decreases with increasing loop order for a fixed $N_{f}$. However, at smaller values of $N_{f}$ there is a discernible deviation amongst the various Padé approximants. A similar phenomenon was apparent for the chiral Ising QED-GNY theory Ihrig et al. (2018a); Zerf et al. (2018) as well as for pure QED Di Pietro and Stamou (2017), whereas Padé approximants for pure GNY models are comparatively better behaved Mihaila et al. (2017); Zerf et al. (2017); Ihrig et al. (2018b). Physically this can be understood from the fact that the disordered phase of pure GNY models (a Dirac semimetal) is adiabatically connected to a system of noninteracting Dirac fermions regardless of the value of $N_{f}$, whereas in QED-GNY models the disordered phase (the ASL) consists of a system of mutually interacting Dirac fermions and gauge fields which becomes increasingly strongly coupled in the infrared for small $N_{f}$ (at least in the sense of the $1/N_{f}$ expansion). In particular, the large variation of the [1/3] Padé approximant for both $\eta_{\phi}$ and $\omega$ in relation to the other approximants was a feature also noticed for the chiral Ising QED-GNY model Zerf et al. (2018); if this approximant is ignored for those exponents a better overall consistency is achieved. Unitarity bounds Ferrara et al. (1974); Mack (1977) in conformal field theory require the scaling dimension $\Delta$ of a Lorentz scalar to obey $\Delta\geq d/2-1$. Since $\Delta_{\phi}=(d-2+\eta_{\phi})/2$ and $\Delta_{\phi^{2}}=d-1/\nu$, this condition imposes $\eta_{\phi}\geq 0$, which is satisfied by all Padé approximants, and $1/\nu\leq d/2+1$. While all the approximants for $1/\nu$ satisfy the latter criterion in $d=3$, some of them give unphysical negative extrapolation values for small values of $N_{f}$.

For fixed fermion flavor number $N_{f}=1$, the behavior of the Padé approximants as functions of $\epsilon$ is shown in Figs. 5-7. As expected, one can observe that all approximants converge to the same result in the $\epsilon\rightarrow 0$ limit. The $\epsilon$ dependence is beneficial to illustrate (1) that there are no singularities in the extrapolation region $0<\epsilon<1$, otherwise the extrapolation would not be well defined, and (2) to show that certain extrapolations violate unitarity bounds before the $\epsilon\rightarrow 1$ limit is reached. In Fig. 5 there is significant variation—as much as two orders of magnitude—among the various approximants at $\epsilon=1$ for $\eta_{\phi}$, although the [1/1] and [2/2] approximants produce similar estimates. (The [2/2] approximant was not plotted in Fig. 2 as it contained singularities in the extrapolation region $0<\epsilon<1$ for certain values of $N_{f}$ in the range considered.) In contrast to the present results, for the case of the pure chiral Heisenberg GNY model, where the gauge field is absent, the Padé approximants for $\eta_{\phi}$ shown in Ref. Zerf et al. (2017) for small $N_{f}$ were in good agreement with each other as well as with QMC and functional RG results. This suggests that gauge fluctuations are a decisive feature in the critical behavior of the scalar field. The stability exponent $\omega$ exhibits less variability than $\eta_{\phi}$, as seen in Fig. 7, but better conformity between approximants was again exhibited in the pure chiral Heisenberg GNY case Zerf et al. (2017). This reiterates the point that for small $N_{f}$ gauge fluctuations tend to destabilize the system. The distinctive behavior of the [1/3] approximant at small $N_{f}$ for $\eta_{\phi}$ and $\omega$, already discussed above, is again illustrated by its unique nonmonotonic dependence on $\epsilon$ at $N_{f}=1$ in Figs. 5 and 7, which was also observed in the chiral Ising QED-GNY case (see Ref. Zerf et al. (2018), as well as unpublished results for $\omega$ vs $\epsilon$). Dropping this approximant reduces the range of extrapolated values for those exponents, but not sufficiently so as to allow us to produce a meaningful quantitative estimate. By contrast, the three approximants [0/2], [0/3], [1/3], while at different loop orders, produce a reasonably consistent estimate of $1/\nu$ for $N_{f}=1$ in the range $\sim$ $0.1$ - $0.3$ (Fig. 6). Interestingly, in the chiral Ising QED-GNY model the same set of three approximants also produced positive values of $1/\nu$ at $N_{f}=1$, and in the same order ($(1/\nu)_{[0/2]}>(1/\nu)_{[1/3]}>(1/\nu)_{[0/3]}$), but spread over a wider range ($\sim$ $0.05$ - $0.7$). Similarly, the [2/1], [2/2], and [3/1] approximants produce unphysical negative values of $1/\nu$ at $\epsilon=1$ in both models.

Turning now to the scaling dimensions $\Delta_{\overline{\Psi}\Psi}$ and $\Delta_{i\overline{\Psi}\Gamma_{5}\Psi}$ of fermion bilinears, corresponding for $N_{f}=1$ to CDW and VBS susceptibility critical exponents, respectively, the Padé approximants are shown for fixed $d=3$ and arbitrary $N_{f}$ in Figs. 8-9, and for $N_{f}=1$ and arbitrary $\epsilon$ in Figs. 10-11. In the large-$N_{\!f}$ limit both scaling dimensions asymptote to $2$, however, for small values of $N_{\!f}$ there is a discernible variation in the approximants as for the exponents considered previously. As we only consider scalar or pseudoscalar bilinears, unitarity bounds require $\Delta_{\overline{\Psi}\mathcal{M}\Psi}\geq 1/2$, which is satisfied by all approximants for $N_{f}\geq 2$, but violated for $N_{f}=1$ by the [1/2] and [0/4] approximants in Fig. 8 (see also Fig. 10 for the [0/4] approximant) and all four-loop approximants in Fig. 9. However, some of those approximants do satisfy the unitarity requirement for $N_{f}=1$ after Borel resummation (see Table 1). There is a higher degree of convergence between approximants at a given loop order for the axial mass bilinear as compared to the normal mass bilinear, as can be seen by comparing Fig. 8 and Fig. 9. Furthermore, the [0/2] (two-loop), [0/3] and [1/2] (three-loop), and [1/3] (four-loop) approximants all give very similar values for the axial mass bilinear dimension at $N_{f}=1$ (Fig. 11), in the range $\sim$ $0.7$ - $0.8$. (The [1/3] approximant was excluded from Fig. 9 for the same reason as that mentioned for $\eta_{\phi}$.)

In Table 1 we collect the numerical extrapolated values for the critical exponents $\eta_{\phi}$, $1/\nu$, and $\omega$, and for the fermion bilinear dimensions $\Delta_{\overline{\Psi}\Psi}$ and $\Delta_{i\overline{\Psi}\Gamma_{5}\Psi}$, as obtained from the $\epsilon$-expansion results using both Padé and Padé-Borel resummation, for the particular case $N_{f}=1$. A high variability in the extrapolation values for $\eta_{\phi}$ is manifest, as already discussed, whereas the variability is significantly less for the other exponents. Supplementing the nonnegative Padé estimates for $1/\nu$ with Padé-Borel resummation extends the range of obtained values somewhat, predicting $1/\nu$ in the range $\sim$ $0.1$ - $0.6$. Similarly, the Padé-Borel values for the fermion bilinear dimensions are systematically higher than the corresponding Padé estimates. Overall, the range of extrapolated values for the exponents, apart from $\eta_{\phi}$, is comparable to that found in the chiral-Ising QED-GNY case. The use of more sophisticated resummation methods Mihaila et al. (2017) may potentially improve those estimates, but typically requires the introduction of a number of free parameters which effectively adds to the uncertainty. Another potentially beneficial approach to improve exponent estimates in three dimensions is the use of two- and four-dimensional perturbation theory Ihrig et al. (2018b); however, at present there is a lack of two-dimensional studies on the chiral-Heisenberg QED3 model.

In Table 2 we present extrapolated values of critical exponents as obtained from the large-$N_{\!f}$-expansion results, again for the particular case $N_{f}=1$. In addition to the critical exponents $\eta_{\phi}$ and $1/\nu$ and the fermion bilinear dimensions considered in Table 1, we also present extrapolated results for the scaling dimension $\Delta_{i\overline{\Psi}\Gamma_{3}\Gamma_{5}\Psi}$ which corresponds to the QAH susceptibility exponent. The leading-order approximants for $1/\nu$ fall within the range predicted by resummation of the $\epsilon$-expansion result. The fermion bilinear results are in good agreement with one another and their uniformity is slightly better than that observed in the chiral-Ising QED3-GNY case Gracey (2018b). The extrapolated values for $\Delta_{\overline{\Psi}\Psi}$ are highly uniform and fall within the broad range of values predicted by resummation of the $\epsilon$ expansion; by contrast, apart from the [0/2] Padé approximant, the large-$N_{f}$ based extrapolated values for $\Delta_{i\overline{\Psi}\Gamma_{5}\Psi}$ are higher than the values predicted by the $\epsilon$-expansion resummation of this quantity. None of the possible approximants for $\eta_{\phi}$ at this order exist, and thus no quantitative statements about this exponent can be made.

V.2 Pure chiral Heisenberg GNY model

We finally turn to the resummation of $\epsilon$-expansion and large-$N_{f}$ expressions for the fermion bilinear scaling dimensions in the pure chiral Heisenberg GNY model. Figures 12 and 13 show the results of Padé extrapolation of the $\epsilon$-expansion expressions for the normal and axial mass operator dimensions, respectively, as a function of $N_{f}$. The corresponding large-$N_{f}$ expressions, Eqs. (139) and (140), are plotted on the same graphs for comparison purposes. The spread of values for $\Delta_{\overline{\Psi}\Psi}$ is appreciable, but comparatively less than for the same quantity in the presence of the gauge field (Fig. 8). This state of affairs was also observed in the chiral Ising QED-GNY model Zerf et al. (2018). Remarkably, the large-$N_{f}$ result and the [0/4] Padé approximant give very similar results for $N_{f}\geq 2$. The spread of values is even smaller, and the convergence of the approximants with increasing loop order better, for $\Delta_{i\overline{\Psi}\Gamma_{5}\Psi}$.

In Figs. 14 and 15 we plot the Padé approximants for $N_{f}=1$ as a function of $\epsilon$, and in Table 3 we present numerical results for the $N_{f}=1$ Padé and Padé-Borel approximants in $d=3$. With the exception of the [0/4] approximant without Borel resummation, the results give a reasonably consistent estimate for $\Delta_{\overline{\Psi}\Psi}$ in the range $\sim$ $2.5$ - $2.7$, and for $\Delta_{i\overline{\Psi}\Gamma_{5}\Psi}$ in the range $\sim$ $1.2$ - $2.0$. For comparison, in Table 4 we display the Padé and Padé-Borel approximants obtained from the large-$N_{f}$ expansion results. While the approximants for the axial mass dimension fall within the range predicted by the $\epsilon$-expansion results, those for the normal mass dimension fall either within the same range (the [0/1] approximant) or give somewhat lower values ($\sim$ $2.1$ - $2.3$).

VI Conclusion

In this paper we provided a comprehensive analysis of the critical properties of the chiral Heisenberg QED3-GNY model, motivated by a recent QMC study Xu et al. (2019) of a $U(1)$ lattice gauge theory with spinful fermions on the 2D square lattice that observed a direct transition between a $U(1)$ deconfined phase, adiabatically connected to the ASL, and an antiferromagnetic Néel-ordered phase. Using the $\epsilon$ expansion below four spacetime dimensions to four-loop order, we showed the existence of a bona fide critical fixed point in the chiral Heisenberg QED3-GNY model for an arbitrary number $N_{f}$ of $SU(2)$ doublets of four-component Dirac fermions. The existence of a critical fixed point for $N_{f}=1$ establishes that the Néel-ASL transition should be continuous, in agreement with the numerical results in Ref. Xu et al. (2019). This fixed point was also shown to exist in the large-$N_{f}$ expansion in fixed $2<d<4$ spacetime dimensions. The $\epsilon$-expansion was used to compute several critical exponents to $\mathcal{O}(\epsilon^{4})$, including the order parameter anomalous dimension $\eta_{\phi}$, the inverse correlation length exponent $1/\nu$, and the stability critical exponent $\omega$. We additionally computed the scaling dimensions of two Dirac fermion bilinears, the normal mass operator $\overline{\Psi}\Psi$ and the axial mass operator $i\overline{\Psi}\Gamma_{5}\Psi$, which for $N_{f}=1$ had the interpretation of CDW and VBS susceptibility exponents at the Néel-ASL transition. These $\epsilon$-expansion results were supplemented by computations in the large-$N_{f}$ expansion. The exponents $\eta_{\phi}$, $\Delta_{\overline{\Psi}\Psi}$, and $\Delta_{i\overline{\Psi}\Gamma_{5}\Psi}$ were computed to $\mathcal{O}(1/N_{f}^{2})$, while $1/\nu$ was computed to $\mathcal{O}(1/N_{f})$. Additionally, the QAH susceptibility exponent, controlled by the scaling dimension of the time-reversal- and parity-odd bilinear $i\overline{\Psi}\Gamma_{3}\Gamma_{5}\Psi$, was computed to $\mathcal{O}(1/N_{f})$. The CDW and VBS exponents were also computed at the critical point of the pure chiral Heisenberg GNY model, which for $N_{f}=1$ describes the semimetal-AF insulator transition in the honeycomb lattice and the $\pi$-flux square lattice. Padé and Padé-Borel resummation techniques were subsequently applied to obtain numerical estimates of all critical exponents in $d=3$ for general $N_{f}$, with a special emphasis on $N_{f}=1$. The critical exponents computed here are in principle accessible to sign-problem-free QMC simulations Xu et al. (2019), and we hope to compare the estimates obtained here to numerical studies of critical properties at the Néel-ASL QCP in the near future. Regarding possible experimental realizations, besides solid-state frustrated magnets we take note of the recent progress in the quantum simulation of lattice gauge theories Zohar et al. (2015). (1+1)D QED was recently simulated experimentally in an ion-trap setup Martinez et al. (2016), and concrete proposals exist to engineer lattice QED in (2+1)D Zohar et al. (2013). The latter would represent a direct physical realization of Hamiltonian (II.1).

Besides its applications to condensed matter physics, the critical fixed point of the chiral Heisenberg QED3-GNY model is of interest as an example of (2+1)D conformal field theory. Recent years have witnessed a resurgence of interest in such theories, due in large part to highly successful numerical implementations of the conformal bootstrap program (for a recent review, see, e.g., Ref. Poland et al. (2019)). These have led to the determination of $d=3$ critical exponents with unprecedented accuracy in various models of interest to both statistical mechanics/condensed matter physics and high-energy physics, such as the 3D Ising El-Showk et al. (2012) and $O(N)$ vector Kos et al. (2014) models, the chiral Ising GNY model Iliesiu et al. (2016), and QED3 Chester and Pufu (2016b); Li (2018). We hope that this work will stimulate the study of chiral Ising/XY/Heisenberg QED3-GNY models with the conformal bootstrap technique, which may lead to precise determinations of critical exponents that could be compared with $\epsilon$- and large-$N_{f}$ expansion results presented in this and previous work.

Acknowledgements.

We gratefully acknowledge I. Affleck, É. Dupuis, S. Giombi, L. Janssen, I. R. Klebanov, X.-Y. Song, K. Wamer, and W. Witczak-Krempa for helpful discussions. RB was supported by the Theoretical Physics Institute at the University of Alberta. The work of JAG was supported by a DFG Mercator Fellowship and he thanks the Mathematical Physics Group at Humboldt University, Berlin for its hospitality where part of this work was carried out. JM was supported by NSERC grant #RGPIN-2014-4608, the CRC Program, CIFAR, and the University of Alberta. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 764850 (SAGEX).

Appendix A Symmetry transformation properties of fermion bilinears

In this Appendix we determine how Dirac fermion bilinears of the form $\overline{\psi}\tau_{i}\sigma_{j}\psi$, that are invariant under $d=3$ Lorentz transformations in the continuum theory, transform under the symmetries of the microscopic Hamiltonian (II.1). (The $\tau_{i}$ matrices act in nodal ($\pm$) space, while the $\sigma_{j}$ matrices act in spin space.) This in turn allows us to associate to each bilinear a microscopic operator with the same symmetry properties, but defined in terms of the original lattice fermions $c_{r\sigma},c_{r\sigma}^{\dagger}$, whose long-distance correlations at the Néel-ASL QCP will be governed by the scaling dimension of the bilinear at the chiral Heisenberg QED3-GNY fixed point.

The microscopic symmetries of interest are the symmetries of the $p4m$ space group of the square lattice, generated by the four-fold rotation $C_{4}$ about a site and the reflection $R_{x}$ about the $yz$ plane and through a site, the unit lattice translations $T_{x}$ and $T_{y}$, and time-reversal symmetry $\mathcal{T}$. The model also has a particle-hole symmetry at half filling, but we will not discuss it further. Choosing the two Bravais lattice vectors for the $\pi$-flux phase in Fig. 1 as ${\bm{a}}_{1}=\hat{{\bm{x}}}-\hat{{\bm{y}}}$ and ${\bm{a}}_{2}=\hat{{\bm{x}}}+\hat{{\bm{y}}}$, and denoting the annihilation operator $c_{{\bm{R}}A\sigma}$ ($c_{{\bm{R}}B\sigma}$) for a fermion of spin $\sigma$ on Bravais lattice site ${\bm{R}}=n_{1}{\bm{a}}_{1}+n_{2}{\bm{a}}_{2}$, $n_{1},n_{2}\in\mathbb{Z}$ and sublattice A (sublattice B) by $c_{(n_{1},n_{2})A\sigma}$ ($c_{(n_{1},n_{2})B\sigma}$), the fermion annihilation operator transforms as

[TABLE]

as well as $c_{(n_{1},n_{2})A\sigma}\rightarrow(i\sigma_{2})_{\sigma\sigma^{\prime}}c_{(n_{1},n_{2})A\sigma^{\prime}}$ and $i\rightarrow-i$ under $\mathcal{T}$, and likewise on the B sublattice.

To determine how the Dirac fermion bilinears $\overline{\psi}\tau_{i}\sigma_{j}\psi$ transform under the microscopic symmetries, one must first determine the PSG Wen (2002) associated with the Hamiltonian (2), i.e., the combinations of transformations (A.1)-(A.8) and $U(1)$ gauge transformations

[TABLE]

that leave the $\pi$-flux Hamiltonian (2) invariant. Writing any of the previously listed symmetry transformations as $\mathcal{S}c_{r\sigma}\mathcal{S}^{-1}=S_{\sigma\sigma^{\prime}}c_{S(r)\sigma^{\prime}}$, where $S(r)$ denotes the transformed coordinate, the associated PSG transformation is

[TABLE]

where $G^{S}_{r}\in U(1)$. By explicit calculation one finds that the gauge transformations in the PSG can be chosen as

[TABLE]

Next one can determine how the PSG acts on the Dirac fermions by following the procedure detailed in Ref. Hermele et al. (2008), i.e., by expanding the Dirac fermion fields in the ${\bm{q}}\rightarrow 0$ limit in terms of the microscopic fermion operators, and using Eq. (A.11). We find

[TABLE]

where the $\mu_{i}$ act on the two-component spinor (Dirac) indices and $K$ denotes complex conjugation. One can check that these are indeed symmetries of the Dirac Hamiltonian (8), provided the momentum ${\bm{p}}$ is appropriately transformed. From Eq. (A.17)-(A.21) we can in turn determine the PSG transformation properties of spin-triplet (Table 5) and spin-singlet (Table 6) Dirac fermion bilinears.