Monte Carlo Study of Compact Quantum Electrodynamics with Fermionic Matter: the Parent State of Quantum Phases

The interplay between lattice gauge theories and fermionic matter accounts for fundamental physical phenomena ranging from the deconfinement of quarks in particle physics to quantum spin liquid with fractionalized anyons and emergent gauge structures in condensed matter physics. However, except for certain limits (for instance large number of flavors of matter fields), analytical methods can provide few concrete results. Here we show that the problem of compact $U(1)$ lattice gauge theory coupled to fermionic matter in $(2+1)$D is possible to access via sign-problem-free quantum Monte Carlo simulations. One can hence map out the phase diagram as a function of fermion flavors and the strength of gauge fluctuations. By increasing the coupling constant of the gauge field, gauge confinement in the form of various spontaneous symmetry breaking phases such as valence bond solid (VBS) and N\'eel antiferromagnet emerge. Deconfined phases with algebraic spin and VBS correlation functions are also observed. Such deconfined phases are an incarnation of exotic states of matter, $i.e.$ the algebraic spin liquid, which is generally viewed as the parent state of various quantum phases. The phase transitions between deconfined and confined phases, as well as that between the different confined phases provide various manifestations of deconfined quantum criticality. In particular, for four flavors, $N_f = 4$, our data suggests a continuous quantum phase transition between the VBS and N\'{e}el order. We also provide preliminary theoretical analysis for these quantum phase transitions.

The interplay between lattice gauge theories and fermionic matter accounts for fundamental physical phenomena ranging from the deconfinement of quarks in particle physics to quantum spin liquid with fractionalized anyons and emergent gauge structures in condensed matter physics. However, except for certain limits (for instance large number of flavors of matter fields), analytical methods can provide few concrete results. Here we show that the problem of compact U(1) lattice gauge theory coupled to fermionic matter in (2 + 1)D is possible to access via sign-problem-free quantum Monte Carlo simulations. One can hence map out the phase diagram as a function of fermion flavors and the strength of gauge fluctuations. By increasing the coupling constant of the gauge field, gauge confinement in the form of various spontaneous symmetry breaking phases such as valence bond solid (VBS) and Néel antiferromagnet emerge. Deconfined phases with algebraic spin and VBS correlation functions are also observed. Such deconfined phases are an incarnation of exotic states of matter, i.e. the algebraic spin liquid, which is generally viewed as the parent state of various quantum phases. The phase transitions between deconfined and confined phases, as well as that between the different confined phases provide various manifestations of deconfined quantum criticality. In particular, for four flavors, N f = 4, our data suggests a continuous quantum phase transition between the VBS and Néel order. We also provide preliminary theoretical analysis for these quantum phase transitions.
In recent years, collective efforts from both condensed matter and high energy physics have started to generate promising * wanderxu@gmail.com † yangqi137@icloud.com ‡ zymeng@iphy.ac.cn outcomes [4,12,[17][18][19][20]. There have been concrete examples, by now, of discrete Z 2 gauge field theories coupled to fermionic matter at (2 + 1)D, deconfined phase with fractionalized fermionic excitations at weak gauge fluctuation, as well as symmetry breaking phase with gapped fermionic excitations at strong gauge fluctuation have been observed [18][19][20]. The apparently continuous transition between deconfined and confined phases is highly non-trivial [20] as it is driven by the condensation of emergent fractionalized excitations and is hence beyond the scope of the Landau-Ginzburg-Wilson paradigm of critical phenomena in which symmetry-breaking is described by a local order parameter. The cQED 3 is the simplest theory to discuss confinement and chiral symmetry breaking [21][22][23]. The pure cQED 3 without matter field is known to be always confining [10,21,24,25]. However, when there is fermionic matter, the gapless fermionic fluctuations may drive the system towards deconfinement. The large N f limit of cQED 3 with fermionic matter is believed to belong to this case [11,13,15,26]. The finite N f case is uncontrolled analytically while numerical methods still face much difficulties [12,17,[27][28][29]. cQED 3 with finite N f flavors of fermionic matter is particularly important to condensed matter physics because these cases actually correspond to the low energy field theory description of many interesting strongly correlated electron systems, and therefore host the potential promise of establishing the new paradigms in condensed matter physics.
Based on these considerations, in this work, we succeeded in performing large-scale quantum Monte Carlo (QMC) simulations on the cQED 3 coupled to N f -flavor of fermions, and eventually map out the phase diagram ( Fig. 1) in the fermion flavor and gauge field fluctuations strength plane. Deconfined phases -U(1) deconfined phase (U1D hereafter) to be more precise -are indeed found in the phase diagram for N f = 6 and 8, and even at N f = 2 and 4 there are very positive signatures of their existence. Various confined phases, in the form of different symmetry breakings, such as antiferromagnetic order (AFM) and valence bond solid (VBS), are also discovered. Interesting quantum phase transitions, between deconfined and confined phases, and between different confined phases [18,[30][31][32][33][34] are revealed as well.
For the sake of smoother narrative, the rest of the paper is organized in the following order. In Sec. II A and II B, we first start with a quantum rotor model coupled to fermions, which can be formulated as cQED 3 coupled to fermionic matter. Then in Sec. II C, we discuss the sign structure of this model, where we find that a pseudo-unitary group can be used to avoid the phase problem at odd N f , and the sign problem at even N f . In Secs. II D and II E, we explain the challenges in the QMC simulation even without sign-problem and provide our solution with a fast update algorithm for simulating gauge fields with continuous symmetries. In Sec. III we discuss the whole phase diagram, and then focus on the physical properties and understanding of the U1D phase, in particular the reason for it being the parent state of various quantum phases, the deconfinement-confinement phase transitions, and VBS to AFM phase transition at N f = 4. Preliminary theoretical analysis of these transitions is also given in Sec. III. Finally, the discussion and conclusions are given in Sec. IV.

A. Rotor model with fermion
The system we are interested in, can be most conveniently formulated as a 2D quantum rotor model coupled to fermions with Hamiltonian whereL i j andθ i j are canonical angular momentum, [L i j , e ±iθ i j ] = ±e ±iθ i j , and its coordinate operator of rotors on each bond b = i, j of a 2D square lattice, as depicted in Fig. 1 (b). The fermion flavor α runs from 1 to N f and the fermions are minimally coupled to the rotor via nearest-neighbor hopping on the square lattice. The flux term with K > 0 favors π-flux in each elementary plaquette , where the magnetic flux of each plaquette is defined as curlθ = b ∈ θ b and the summation overθ b has been taken in either clockwise or anticlockwise orientation around an elementary plaquette.
For the Monte Carlo simulations, it is convenient to work in a representation whereθ i, j is diagonal. That is, omitting the bond index,θ|φ = φ|φ with φ ∈ [0, 2π]. In this representation,L = −i ∂ ∂φ with eigenvectorsL φ|l ≡Le iφl = le iφl and l ∈ Z. With these definitions, the resolution of unity reads: ∫ 2π 0 |φ φ| = 1 2π l |l l | =1. To formulate the path integral, we have to estimate the matrix element: φ |e −∆τJ N fL 2 /8 |φ . To this end we insert resolution of the unit operator, and use the Poisson summation formula to obtain: The last equality corresponds to the Villain approximation.
With the above, the Hamiltonian in Eq. (1) can be formulated in a coherent-state path integral with action S = S F + S φ = ∫ β 0 dτ(L F + L φ ) and the Lagrangian for fermion and gauge field parts are respectively, where µ will be set to zero for half-filled case. β = 1 T is the inverse temperature. The model in Eq. (1) has now been explicitly formulated as (unconstrained) cQED 3 coupled to fermionic matter [11,13,17]. We will now consider the symmetries of the model that will allow to prove that the constraint is dynamically imposed in the low temperature limit.

B. Symmetries and limiting cases
Our model, see Eq. (1), has global and local symmetries. It enjoys a manifest global SU(N f ) spin symmetry as well as a particle hole-symmetry: In the above, z is a complex number that makes it clear that the particle-hole symmetry is anti-unitary, and (−1) i takes the value 1 (−1) on sublattice A (B). The local U(1) gauge transformation is an invariant. The generator of this local symmetry corresponds to a local conserved charge (Gauss law) with [Q i ,Ĥ] = 0. In our simulations we sample over allQ i sector, such that our Hamiltonian corresponds to an unconstrained gauge theory. As a consequence, correlation functions of gauge dependent quantities such as the single particle operator are local in space but not it time: ĉ iα (τ)ĉ jα = δ i, j .
Below we argue that the Gauss law constraint is dynamically imposed in the zero temperature limit. At J = ∞,L i j vanishes and charges are completely localized since hopping on a given bond involves excitations of the rotor mode. In this limit charge configurations, corresponding to specific values ofQ i , are degenerate and at any finite value of J the degeneracy will be lifted. The dynamical generation of the termĤ accounts for the lifting of this degeneracy. Note that since P ,Q i = 0 terms containing products of odd numbers of Q i 's are forbidden. The above equation defines a classical model with a finite temperature Kosterlitz-Thouless transition. At zero temperature theQ i are frozen in a given pattern, and the Gauss law is imposed. For J → ∞ our model maps onto an SU(N f ) quantum antiferromagnetic. We again start from the J = ∞ degenerate case, and consider t in second order degenerate perturbation theory. As mentioned above, hopping of a fermion with flavor index α from site i to nearest neighbor site j leaves the rotor in an excited state, associated with energy cost J. The only way to remove this excitation is for a fermion with flavor index α to hop back from site j to site i. These processes are encoded in the SU(N f ) Heisenberg model: i,αĉj,α . In our simulations we have on average N f /2 fermions per site such that the representation of the SU(N f ) group corresponds to the antisymmetric self-adjoint representation (i.e. Young tableau corresponding to a column of N f /2 boxes). On the square lattice and for even values of N f where the negative sign problem is absent, this model has been considered in former auxiliary field QMC simulations [12]. At N f = 2 one finds an antiferromagnetic state and at and beyond N f = 6 a VBS state. In the large-N f limit we recover the Marston-Affleck [2] saddle point accounting for dimerization. At N f = 4 and in the absence of charge fluctuations, Ref. [12] finds no compelling evidence of VBS and AFM orders when considering lattices up to 24 × 24. On the other hand, simulations of the corresponding SU(4) Hubbard model [35] are consistent with an AFM state in the large U limit albeit with decreasing value of the order parameter as a function of U. In our simulations charge fluctuations are present and the phase diagram is consistent with AFM order in the large J limit.
At J = 0,θ i j becomes a classical variable. Even in the absence of the flux term, the coupling to the fermions favors, according to Lieb's theorem [36], a π-flux per plaquette, with associated dynamically generated Dirac dispersion relation of theĉ-fermions. The fate of this state at low values of N f and when gauge fluctuations are accounted for is one of the central aims of our research.

C. Absence of sign problem for even N f
To simulate above model with quantum Monte Carlo method, we start with the partition function As the action of gauge field part S φ = ∫ β 0 dτ L φ with L φ shown in Eq. (3), is always real (thus its exponential is always positive), the sign structure of the Monte Carlo weight will only come from the trace over fermions. To trace out fermion, we first discretize the imaginary time τ = z∆τ (z = 1, 2, · · · , L τ ) where L τ is the total number of time slices (L τ ∆τ = β), then performing the fermion trace, we have with S F = ∫ β 0 dτ L F and L F shown in Eq. (3), and after the discretization of β, B z = e −∆τV z with V z the coupling matrix for each fermion flavor (we have N f in total, therefore there is power N f in Eq. (10)), which only has elements connecting sites between different sublattice of the square lattice. We recognize that such kinds of B z matrices form a pseudo-unitary group SU(n, n) where 2n is dimension of B z (total number of sites). As proved in Appendix A, ∀D ∈ SU(n, n), det(1 + D) ∈ holds. Therefore the fermion weight will always be real for all integer N f and, most importantly, be semipositive-definite for all even value of N f such that QMC simulations can be performed. Although the current paper only focus on the compact U(1) gauge fields, the SU(n, n) group actually allows one to add extra non-Hermitian terms to the model, such as a staggered imaginary chemical potential, that are also sign problem free for even values of N f . It will be very interesting to study the non-Hermitian models and their properties, in the presence of gauge field fluctuations and electronic interactions in future investigations.

D. Difficulties of the QMC simulation
Although there is no sign-problem for even N f , the simulation of the Eq. (1) is by no means simple. Earlier attempts in the high energy community have been devoted to simulate similar models by means of hybird Monte Carlo [4,17,29,37] with (dynamical) mass term. Mass terms are essential to avoid divergences when calculating forces in the realm of Hamiltonian [38] and Langevin [39] dynamics. Mass terms however introduce biases the severity of which have to be a posteriori clarified.
On the other hand, in the condensed matter community, the determinantal QMC (DQMC) is more popular and it usually uses local updates and the mass terms are not essential here [40][41][42][43]. However, as far as we are aware of, there exists no DQMC simulation on cQED 3 coupled to fermionic matter, and our work hence serves as a first attempt. As explained below, to be able to simulate the model in Eq. (1), there are several obstacles one needs to overcome.
The most obvious obstacle is about the computational complexity. For general models the computation complexity of DQMC for one sweep is O(βN 4 ) (here N is total number of sites), for models where fast update is applicable the complexity can be reduced to O(βN 3 ). The most common example is Hubbard model with onsite interaction, when we flip a single auxiliary field at time slice z and site i, it only changes one diagonal value in B z matrix, then the new equal time Green's function G (τ, τ) can be calculated as where G is the Green's function at previous step. As ∆ ii is the only non-zero element in N × N matrix ∆, ∆(1 − G) can be written as the outer product of two vectors, then Sherman-Morrison formula can be used to reduce the complexity for calculating Eq. [40,41]. Such update scheme is referred as fast update. For one sweep over all auxiliary fields (scales as O(βN)), the scaling will be O(βN 3 ) instead of O(βN 4 ). For models with off-site interaction, for example in our case due to the coupling between fermion and gauge field, usually we need to make a further Suzuki-Trotter decomposition over all bonds (assume bond b connects site i and j, total number of bonds is N b ), and B z matrix will be written as where is the complex function of auxiliary field φ i j in polar form. It is obvious that only when θ(φ i j ) = 0 or π, B z,b (i, j) will be diagonalized by an auxiliary field independent unitary transformation, and then the aforementioned fast update scheme applies. There are many known models belong to this case, such as model with Heisenberg type interaction [44][45][46], models with Z 2 (bosonic) gauge field coupled to fermionic matter [18,19,47,48], etc. Unfortunately, our model in Eq. (1) involves U(1) gauge fields, the auxiliary field φ i j is therefore a continuous variable and θ(φ i j ) = φ i j , thus our model does not belong to the cases discussed above and naively the fast update cannot be applied.

E. Fast update algorithm designed for U(1) gauge fields
There is indeed an alternative way to to design a fast update algorithm for the model in Eq. (1). It is based on the Woodbury matrix identity that effectively generalizes the Sherman-Morrison formula to higher rank matrices. We recognize that the new B z,b after a single update can directly be factorized and other elements of the N × N matrix ∆ are zero. With such kind of structure and note that Eq. (11) still holds, ∆(1−G) can be written as product of two matrices with dimension N × 2 and 2 × N. Therefore, we can use the generalized version of Sherman-Morrison formula (Woodbury matrix identity) to calculate Eq. (11), which also has complexity O(N 2 ). With such special designed fast update, we are now ready to simulate cQED 3 coupled to fermionic matter without artificial mass term and still enjoy the O(βN 3 ) computational complexity.

A. Physical observables
Our model and major results are schematically summarized in Fig. 1 (a) and (b), respectively, but before starting the discussion of QMC results, we first introduce the QMC observables that are used to characterize symmetric and symmetry breaking phases and their phase transitions. Since physical observables are Hermitian, we constructed and measured various gauge invariant structure factors, including spin χ S (k) and dimer χ D (k) structure factors. They are defined as where the spin operator S α is defined as dimer along the nearest-neighbor bond inx direction.
From these structure factors, one can further construct dimensionless quantities -the correlation ratio [49] -to determine the precise position of phase transitions. If one would like to detect the transition towards antiferromagnetic long range order, the antiferromagnetic correlation ratio is where X = (π, π) is the order wavevector for AFM on the square lattice and δq = ( 2π L , 0) is the smallest momentum away from X. In the same vein, we define correlation ratio for the valence bond solid (VBS) order from dimer structure factor where M = (π, 0) is the order wavevector for VBS. Other quantities, such as the energy density and various correlation functions (spin, dimer) in real space, are also measured in the QMC simulation.

B. Phase diagram
Now we can discuss the results from QMC simulation of Eq. (1). Starting with the final phase diagram that schematically summarized all the data, as shown in Fig. 1 (a), the phase diagram is spanned along the axes of N f and J. We set K = t = 1 as the energy unit and choose the gauge fluctuation strength J as the tuning parameter to study different N f cases.
For each N f there are different phases and phase transitions, but there are similarities for all N f investigated, that is, at small J, U(1) deconfined phases (U1D in Fig. 1 (a)) are universally present in the phase diagram. This finding is highly non-trivial, as explained in the introduction (Sec. I), from both high energy physics and condensed matter physics communities, the existence of such a deconfined phase in cQED 3 is still under debate, due to the lack of controlled calculation at finite and small N f [4][5][6][7][8][9][10][11][12][13][14][15], and our finding presented here provide the first set of concrete evidence to support the existence of this phase.
Moreover, as will be further elucidated in later sections, such a deconfined phase is expected to be the algebraic spin liquid [6-9, 13, 14, 50], in which many competing orders, such as antiferromagnetic order, valence bond solid order, charge density wave and superconductivity coexist in a critical manner, and share the same power-law decay in the correlation functions due to the U(1) gauge deconfinement and the subsequential conformally invariant, interacting fixed point [13,14]. Starting from the algebraic spin liquid phase, one could easily apply various perturbations and drive the system into various symmetry-breaking phases. Therefore, the algebraic spin liquid phase, i.e. U1D discovered in this work, actually serves as the original state of many interesting quantum phases, hence dubbed parent state of quantum phases. The discovery of such a deconfined phase, is the most important result of this work.
As J increases, the system goes through deconfine-confine phase transitions to various symmetry-breaking phases. In the case of N f = 2, the symmetry-breaking phase is the AFM (Néel) phase, whereas in the case of N f = 4, the symmetrybreaking phases are VBS and AFM phases. And further increases N f , the symmetry-breaking phases are VBS solely. According to Ref. [20], the deconfine-confine phase transition could be a version of the deconfined quantum critical point with emergent continuous symmetry, and besides the QMC data, we will also provide a preliminary field theoretical description of this transition. The transition from VBS to AFM phases inside the confined regime at N f = 4, if it is indeed continuous as our data suggest, is also a deconfined quantum critical point [18, 30-34, 51, 52], whose theory will be explained later.
Overall, the presence of the U1D deconfined phase, and the phase transitions between deconfined to confined phases and within the confined phases, are all intriguing and show the rich physics behind the simple model of Eq. (1).
Below, we discuss some exotic aspects of the phases and phase transitions we obtain.

U1D phase
The U1D phase found in the small J region is expected to have the same power law decay for both the dimer correlation and spin correlation [13]. In fact, according to the large-N f perturbative renormalization group calculation, these correlation functions decay as [7,13,14,50]. Notice that in our case, N f is the number of fermion flavors on the lattice, while in Refs. [7,13,14,50], N f is the number of two-component Dirac fermions, which is twice of our N f due to momentum valley degeneracy. We now compare this theoretical expectation to our numerical simulation results. For N f = 8, which is the largest N f simulated in this paper, the correlations of dimer and spin indeed share the same power-law decay. At N f = 8 a study of the VBS correlation ratio presented in Appendix B shows that J c = 2.5(3) (see Fig. 12). Fig. 2 plots the spin-spin and dimerdimer correlations at J = 2.3. As expected for the algebraic spin liquid phase, both quantities follow the same powerlaw and from the QMC simulation we estimate the exponent to be 3.0(6). It is remarkable to see that this result is comparable with large N f value.
For smaller values of N f the simulations in the U1D phase become progressively more demanding due to long autocorrelation times. As argued in Appendix B we are nevertheless able to estimate J c and we obtain J c = 1.6(3) for N f = 6; J c = 0.44 (29) for N f = 4; J c = 1.43 (46) for N f = 2. A noticeable feature of the data is the increase of J c from N f = 4 to N f = 6 and 8. N f = 2 stands at odds with other flavor values in the sense that the instability of the U1D phase is towards an AFM rather than toward the VBS at N f ≥ 4.
Next, and as argued in Appendix B, we can estimate the power law of spin-spin and dimer-dimer correlations in the UID phase as a function of N f . Due to long autocorrelation times it is hard to carry out simulations deep in the U1D phase and the results presented here, especially for the small N f , are taken close to the finite size value of J c . This quantity corresponds to the crossing point of the correlation ratio between adjacent system sizes. Our results are consistent with identical power laws for dimer and spin correlation functions for each value of N f . Fig. 3 presents a summary plot of the power law we obtain at N f = 8, 6, 4, and 2. It is remarkable to see that our data progressively approaches the aforementioned 1/N f perturbative expression.

Confinement transition
We find a U1D to AFM phase transition for the N f = 2 case and a U1D to VBS phase transition for the N f = 4, 6 and 8 cases. These phase transitions should belong to the QED 3 -Gross-Neveu-O(3) or XY transitions [53], depending upon the order parameters in the confined phases. For example, the transition between U1D to VBS phase can be described by the following action: vector in which the VBS order parameter is embedded. µ = (µ x , µ y ) are two 2N f × 2N f matrices in the fermion flavor space.ψ µ x ψ,ψ µ y ψ are two fermion mass operators that correspond to the VBS in x and y directions respectively. When r > 0, φ is gapped out, and the system is in the U1D phase due to the screening of massless fermions to the gauge field. When r < 0, φ condenses, the fermions are gapped out, then the compact gauge field is in the confined phase. In the U1D phase, the VBS and antiferromagnetic order parameters should have the same scaling dimension, due to the enlarged SU(2N f ) symmetry in the low energy field theory, consistent with our data in Fig. 2 and Figs. 7,9 and 11 in Appendix B; while at the critical point r = 0, these two order parameters still have power-law correlation functions, but with different scaling dimensions, due to the loss of the SU(2N f ) symmetry in the infrared.

AFM-VBS transition at
The situation at N f = 4 is even more interesting. As we further increase J, we observe another quantum phase transition from VBS to AFM phase. This transition is consistently revealed in three steps in Fig. 4. Fig. 4 (a) shows the r AFM correlation ratio and clearly there is a crossing point signifying the establishment of the AFM long range order. The inset shows the 1/L extrapolation of the crossing point and gives rise to J c2 = 18(3). Fig. 4 (b) is the r VBS correlation ratio and the crossing point in it signifies the vanishing of the VBS order. Inset of Fig. 4 (b) gives rise to J c2 = 19(5), consistent with the onset of the AFM order in Fig. 4 (a).
The transition from VBS to AFM deserves more attention, apparently the data in Fig. 4 suggest a continuous transition, and if it were the case, there is then the possibility that the critical point will acquire larger symmetry group than that in the model in Eq. (1). In the case of Z 2 gauge field coupled to fermion, as shown in Refs. [18,20], two similar situations with emergent continuous symmetry are also investigated. In the first case, it is at N f = 3 that a continuous VBS to AFM phase transition occurs [18], and in the second case, it is the deconfine-confine phase transition itself at N f = 2 [20]. Our phase transition at J c2 is closer to the former. To further understand the nature of transition from VBS to AFM, we plot the ratio of AFM structure factor and VBS structure factor, The results are shown in Fig. 4 (c), indeed there is a crossing point in the r AFM/VBS , and the 1/L extrapolation in the inset of   Fig. 4 (a) and r VBS in Fig. 4 (b). This transition is similar to the AFM-to-VBS transition in the SU(4) J-Q model [54], where a continous transition is observed. The transition in that case can be described by a NCCP N −1 description with N = 4 [55], and it is shown that the monopoles are irrelevant at this fix point [56] and therefore a deconfined quantum critical point [30] is realized. However, in Ref. [55], on sublattice A of the lattice there is a fundamental representation of SU(4), while on sublattice B there is a anti-fundamental representation, which is different from our case. In our case, with N f = 4 there is a self-conjugate representation of the SU(4) group on every site, thus the field theory for the AFM-to-VBS transition is different from Ref. [55]. According to Ref. [57], the AFM Néel order in this case has the following grassmanian ground state manifold M: To describe this antiferromagnetic state, one can either introduce N f = 4 flavors of fermionic spinons with halffilling, or introduce two color species of bosonic spinons z α,a (α = 1, · · · 4, a = 1, 2), and couple them to a U(2) gauge field. The U(2) gauge constraint will guarantee that on every site there are fixed number of spinons, and the color space is fully antisymmetric, thus on every site the SU(4) spin is automatically in an antisymmetric self-conjugate representation. Then the field theory for the Néel-VBS transition is where a µ and a l µ are gauge fields corresponding to the U(1) and SU(2) subgroups of U (2). Note that these gauge fields are "emergent" gauge fields, which are different from the explicit gauge field in our original simulation.
When r < 0, z α,a condenses, and leads to the antiferromagnetic state with ground state manifold U(4) U(2)×U (2) . When r > 0, z α,a is gapped out, and the gauge fields will be confined. Here we assume that the U(1) compact gauge field a µ still has the quadru-monopole proliferation, which leads to the VBS phase like the original deconfined QCP theory for the SU(2) spins [30].
One can also describe the SU(4) Néel and VBS order parameters on an equal footing, and capture the "intertwinement" between the two order parameters with a topological term. To do this, we need to embed both the Néel and VBS order parameters into a larger manifold. One way to parameterize the Néel order manifold is N = U † ΩU, where Ω is a 4 × 4 diagonal matrix Ω = diag(1 2×2 , −1 2×2 ), and U is a SU(4) matrix. N is a 4 × 4 Hermitian matrix with constraint N 2 = 1. Now we introduce a larger 8 × 8 matrix P which includes both the Néel and VBS order parameters: where (V x , V y ) is a two component order parameter for the VBS phase, and V 2 x + V 2 y = 1. The order parameter P unifies the SU(4) Néel and VBS order parameters, just like the O(5) vector order parameter introduced in Ref. [58].
The topological term that captures the "intertwinement" between Néel and VBS order parameters is a Wess-Zumino-Witten term: (24) Using the same technique in Ref. [59], one can show that at the vortex core of the VBS order parameter, there is a spinon with self-conjugate representation, which is consistent with intuition. In fact, the O(5) WZW term introduce in Ref. [58] can be written in the same form as Eq. 24, as long as we replace P in Eq. (24) by a 4 × 4 Hermitian matrix order parameter P = n · Γ, where Γ are five Gamma matrices, and n is the O(5) vector introduced in Ref. [58].
This topological term can be viewed as the low-energy effective field theory of the π−flux state of the SU(4) antiferromagnet, which again is described by a QED 3 with eight flavors of Dirac fermions [13], but again the gauge field of this QED 3 is an emergent gauge field which is different from the gauge field introduced in the original model that we simulate. The WZW term Eq. (24) can be derived by coupling the 8 × 8 matrix order parameter P to the eight flavors of Dirac fermions of the π−flux state, and integrate out the fermions following the procedure of Ref. [60].
Our data, the crossing of r AFM/VBS in Fig. 4 (c), suggest that the AFM and VBS order parameters have the same scaling dimension at this transition, which is consistent with the emergent SU(8) symmetry of the π−flux state of the SU(4) antiferromagnet. The large SU(8) symmetry, if indeed exists at the AFM-VBS transition, will ensure that many other order parameters have the same scaling dimension as AFM and VBS order parameters [13]. These order parameters will also have similar fractionalization dynamical signatures in their spectral functions as AFM and VBS order parameters.

IV. CONCLUSIONS
Using large scale DQMC, we investigated the compact U(1) gauge field theory coupled to Dirac fermion matter fields in (2 + 1)D and variable flavor number N f : i.e. cQED 3 . With our simulations we mapped out the entire ground state phase diagram in the flavor N f and gauge field fluctuation J strength plane. Our results are summarized in Fig. 1(a).
Most importantly, signatures supporting stable U(1) deconfined (U1D) phases were discovered at N f = 8 and 6, and evidence of the U1D phase at N f = 4 and 2 were also found. The properties of the deconfined phase are consistent with the proposal of algebraic spin liquid, in which, various competing orders (AFM order and VBS order, for example) all have algebraic correlation with identical power-laws in real-space. The decay power is found to quantitatively converge to the large-N f predictions [7,13,14,50].
The transition between the deconfined and confined phases at various N f were determined using the RG invariant correlation ratios. At N f = 2 the transition occurs between the U1D and AFM phases. Since the AFM corresponds to O(3) symmetry breaking, the critical theory should be described by the QED 3 -Gross-Neveu O(3) universality class. In contrast, at larger values of N f the ordered phases corresponds to VBS. The dynamical generation of the two VBS mass terms is described by the QED 3 -Gross-Neveu O(2) universality class. As far as we know, the QED 3 -Gross-Neveu O(2) or O(3) transition have not been investigated numerically before in an unbiased simulation. It is certainly worthwhile to further carefully study the critical properties of these transitions via QMC simulations and compare with future analytical calculations. In particular, the QED 3 -Gross-Neveu Ising transition has been investigated recently with perturbative RG calculations [53,61].
Aside from the QED 3 -Gross-Neveu transitions, we have found evidence for a direct and continuous transition between the AFM and VBS states in the confined region of the phase diagram at N f = 4. Since we have on average two fermions per site, we should consider the antisymmetric self-conjugate representation of the SU(4) group. We have presented various theoretical descriptions of this putative deconfined quantum phase transition in terms of multi-flavored spinons coupled to emergent U(1) and SU(2) gauge fields. We also discussed the effective low energy field theory with a topological term that captures the intertwinement between the magnetic and VBS orders, and its connection to the π−flux state of the SU(4) spin system discussed in Ref. [13]. Numerical support for emergent symmetries was provided.
Finally, the confining phase transitions in our model, as well as the possible deconfined quantum critical point at N f = 4, will have distinct and very interesting dynamical properties in their spectral functions, that can be further explored in QMC simulations. Such calculations provide experimentally accessible signatures of exotic states of matter where emergent gauge fields, fractionalized excitations, can be traced. Similar attempts have recently been applied to the deconfined quantum critical point in pure spin model [52], emergent Z 2 spin liquid at (2+1)D [62] and U(1) spin liquid at (3+1)D [63] and the Z 2 counterpart of our model [18]. In the present cQED 3 model, dynamical measurements in the QMC simulation plus state-ofart analytical continuation [62,64,65] can help to reveal more fundamental physical understanding of these exotic quantum phase transitions. As we discussed in the main text, the fermion determinant for one flavor is det I n+m + L τ z=1 B z , where n and m are the numbers of sites in the two sublattices, and I n+m denotes the (n + m)-dimensional identity matrix. B z = e h z , where h z has the following structure, and T z is the hopping matrix between different sublattices. B z matrices satisfy (1) B † z ηB z = η, where η = diag(I n , −I m ), and (2) detB z = 1, thus their products generate the pseudounitary group SU(n, m).
Denote the eigenvalues of D by λ i , 1 ≤ i ≤ n + m, then det(I n+m + D) = i (1 + λ i ). We then treat the eigenvalues on the unit circle and those not on the unit circle separately. For those not on the unit circle, (A2) For those on the unit circle, denoting Therefore, we find hence det(I n+m + D) ∈ R. This implies that for an even number of flavors fermions, the model is free of sign problem for any hopping matrices T z . For instance, this is true for both abelian and nonabelian gauge fields.
It is easy to find an example such that det(I + D) < 0. Therefore, the absence of sign problem does not hold for an odd number of fermions in general. However, it does hold for models with fermions coupled to Z 2 gauge fields [18][19][20].
Appendix B: Determination of the U1D phase boundary First we focus on the case N f = 2. In the static limit (J → 0), the gauge fields are frozen into a π-flux pattern per plaquette [36,66] 5 shows the AFM correlation ratio r AFM as a function of J for different system sizes L = 8, 10, 12 and 14. The inset of Fig. 5 traces crossing points of pairs of adjacent system sizes, and a power-law extrapolation in 1/L gives rise to an estimate of the confinement transition in the therodynamic limit: J c = 1.43 (46). Since the correlation ratio shows no abrupt features, the transition from the deconfined phase to the AFM is more likely to be continuous. This is consistent with the flux energy per plaquette data presented in Fig. 13 in Appendix C. For the same J values as considered in Fig. 5 the flux energy per plaquette does not develop a sharp change of slopes. As soon as the matter fields bind to form the particle-hole condensate, we expect monopoles to proliferate and to generate a photon mass. The photon mass can be measured by the correlation of flux quantity θ(τ) [67], which is defined as θ(τ) = sin (curlφ(τ)). The photon mass m is related to the correlation of θ(τ) by C(τ) = θ(τ 0 )θ(τ 0 + τ) ∼ exp(−mτ). Fig. 6 depicts the semi-log and log-log plot of C(τ) for L = 12. On this finite sized system we observe a growth of the photon mass as J increases.
To further understand the properties of the deconfined phase, we measure the real-space correlation functions at J = 1.75. This value of J is larger than our estimate of J c in the thermodynamic limit. However, for finite size systems L = 12 and 14, J = 1.75 is on the U1D phase side of finite size crossing points. As shown in Fig. 7, the spin-spin correlation seems to develop long-range order (decay to a constant value) but the the short range part shows a power-law close to r −1.5 . Interestingly the dimer-dimer correlation function also decay with a power-law close to r −1.5 (dashed line in Fig. 7). This result sheds light on the property of the deconfined phase, which is proposed in Refs. [7,13] to correspond to the algebraic spin liquid. It has the unique property that as a deconfined state emerging from competing orders, the correlation functions of these competing orders, such as spin-spin, dimer-dimer and bond-bond, have the the same power-law decay. If the data in Fig. 7 were deep inside the confined phase, the decay of spinspin and dimer-dimer correlations will be very different. For example, in the Néel phase, spin-spin will decay to a constant value and dimer-dimer correlation will decay exponentially. Therefore, our data in Fig. 7 provide supporting evidence of the algebraic spin liquid behavior of the U1D in Fig. 1 at Next we turn to the N f = 4 case where we also observe a U1D phase at small J. As shown in Fig. 8, we can follow the crossing points of the correlation ratio of the VBS order parameter for different system sizes, so as to extract (see inset) J c1 = 0.44 (29). The data is consistent with a continuous phase transition from the deconfined phase to the VBS phase. Furthermore, the flux energy per plaquette also supports a continuous transition.   Fig. 7, the data are more scattered, but it is still possible to assign a power law decay of the correlations following a ∼ r −2 i j law both for the spin-spin and dimer-dimer correlations. This power law is faster than at N f = 2, and is hence consistent with the large-N f prediction. In particular in the limit N f → ∞ both spin-spin and dimer-dimer correlations are expected to decay as ∼ r −4 i j . Taking into account 1/N f corrections gives  Besides the power-law decay of various competing correlation functions, the situation at N f = 4 is even more interesting than that at N f = 2. As we further increase J, we observe another phase transition from VBS to AFM. This phase transition is discussed in detail in Sec. III B 3 of the main text.
c. N f = 6 and N f = 8 In this section we discuss the results at N f = 6 and N f = 8. As shown in Fig. 10 and Fig. 12, and the corresponding insets, we estimate J c = 1.6(3) for N f = 6 and J c = 2.5 (3) for N f = 8. Again, the data are consistent with a continuous transition between the deconfined UID and confined VBS phases. Correspondingly, the flux energy per plaquette behaves as smooth function across the critical point. More interestingly, we plot the spin-spin and dimer-dimer correlation function in real-space for N f = 6 at J = 1.8 in Fig. 11 and for N f = 8 at J = 2.3 in Fig. 2 of the main text, respectively. In Fig. 11, the spin-spin and dimer-dimer correlation functions show the same power-law decay with exponent 2.8. In Fig. 2, both correlation functions decay with exponent 3.0. This increase in the exponent is consistent with the perturbative RG results. On the whole our data provides concrete evidence that the deconfined phase in our model at various values of N f belong to the algebraic spin liquid [7,13,14,50]. One can foresee that with a further increase of N f , we will reach the expected power-law behavior ∼ r −4 . To characterize the continuous nature of confined and deconfined phase transition, we also measured the flux energy per plaquette 1 L 2 cos curlθ . Fig. 13 depicts our result at N f = 2. For other N f 's, the flux energy per plaquette has a similar continuous behavior.

Appendix D: Dynamically generated constraint
As mentioned in the main text, our model corresponds to an unconstrained gauge theory. As such, the Gauss law will be dynamically imposed andQ i , defined in Eq. 6, should converge to constant value in the zero temperature limit. We studied the uniform structure factor ofQ i by calculating We find that the uniform structure factor ofQ i defined above extrapolates to zero in thermodynamic limit as showed in Fig. 14 for N f = 2. Other N f cases show similar behavior. FIG. 14. Uniform structure factor ofQ i (defined in Eq. D1) for N f = 2. Note that β = 4L here. For all Js, this uniform structure factor extrapolated to zero in thermodynamic limit. Thus, the constrain of Q i = 0 is dynamically enforced.