A simple fermionic model of deconfined phases and phase transitions

Using Quantum Monte Carlo simulations, we study a series of models of fermions coupled to quantum Ising spins on a square lattice with $N$ flavors of fermions per site for $N=1,2$ and $3$. The models have an extensive number of conserved quantities but are not integrable, and have rather rich phase diagrams consisting of several exotic phases and phase transitions that lie beyond Landau-Ginzburg paradigm. In particular, one of the prominent phase for $N>1$ corresponds to $2N$ gapless Dirac fermions coupled to an emergent $\mathbb{Z}_2$ gauge field in its deconfined phase. However, unlike a conventional $\mathbb{Z}_2$ gauge theory, we do not impose the `Gauss's Law' by hand and instead, it emerges due to spontaneous symmetry breaking. Correspondingly, unlike a conventional $\mathbb{Z}_2$ gauge theory in two spatial dimensions, our models have a finite temperature phase transition associated with the melting of the order parameter that dynamically imposes the Gauss's law constraint at zero temperature. By tuning a parameter, the deconfined phase undergoes a transition into a short range entangled phase, which corresponds to N\'eel/Superconductor for $N=2$ and a Valence Bond Solid for $N=3$. Furthermore, for $N=3$, the Valence Bond Solid further undergoes a transition to a N\'eel phase consistent with the deconfined quantum critical phenomenon studied earlier in the context of quantum magnets.

Ground states of strongly interacting electronic systems can exhibit an extremely rich variety of phases. One way to organize our understanding is to classify the zero temperature phases by their entanglement structure at long distances and low energies [1][2][3][4][5][6][7][8][9]. Gapped phases which possess a local order parameter are characterized by short-range entanglement, that is, the reduced density matrix of a large subsystem A can be understood simply by patching together density matrices of smaller subsystems A i whose union is A [10]. This is no longer true for gapless phases such as Fermi liquids [11][12][13], or gapped topological phases such as a fractional quantum Hall liquid, and such phases are thus said to possess 'longrange entanglement' [7,8,14].
Experience as well as heuristic arguments suggest that Hamiltonians whose ground states posess long range entanglement are relatively difficult to simulate on a classical computer. For example, even a phase as ubiquitous and as well understood as a Fermi liquid is rather hard to simulate numerically because fermions at finite density with repulsive interactions tend to have an intricate sign structure in their wavefunctions leading to the infamous Monte Carlo fermion sign problem [15][16][17]. Similarly, fractional quantum Hall phases again possess a non-trivial sign structure to their wavefunctions [18], and therefore so far are amenable only via techniques such as Exact Diagonalization [19,20] and Density Matrix RG [21], which are restricted to only small twodimensional systems due to exponential scaling of numerical cost with system size. However, there do exist long-range entangled phases that do not suffer from Monte Carlo sign problem. Two prominent examples are: (1) Interacting Dirac fermions with an even number of flavors [22] (2) A gapped Z 2 topological ordered system such as a Toric code Hamiltonian [14] in a magnetic field [23]. The absence of sign problem for the former is related to the positive fermion determinant, while in the latter case, the Hamiltonian has nonpositive off-diagonal elements in a local basis, allowing one to sample the corresponding Boltzmann weight efficiently [24]. In this paper, we study a sign problem free model that has a ground state with features of both of the aforementioned long range entangled phases together, namely, Dirac fermions coupled to a fluctuating Z 2 gauge field. By tuning parameters in the Hamiltonian, we also study competition with symmetry breaking short-range entangled phases which leads to novel strongly interacting quantum critical phenomena. We observe a Z 2 Dirac deconfined phase (Z 2 D), a spin-density wave (SDW) phase (≡ Néel phase) (or a superconductor (SC), depending on the pattern of particle-hole symmetry breaking), a charge density wave (CDW) phase as well as a valence bond solid (VBS). For N = 1, we do not find evidence for a Z 2 D phase beyond h = 0 consistent with the arguments in the main text. The phase transitions from the Z 2 D to SDW/SC (N = 2) and VBS (N = 3) are seemingly continuous. At N = 3 we observe a deconfined quantum critical point (DCQP) between the VBS and SDW phases.
Our Hamiltonian (Eq.1) is partially motivated by a recent study of nematic instability of a Fermi surface [25] and velocity fluctuations in Dirac metals [26]. It consists of N species of fermions on the vertices of a square lattice interacting with quantum Ising spins that live on the links of the same square lattice. However, in contrast to the Hamiltonian in Ref. [25], our Hamiltonian has a local Z 2 invariance, which leads to L x L y conserved operators on a L x × L y lattice. Furthermore, we restrict ourselves to half-filling and thus at low energies, instead of a Fermi surface, we either obtain either no fermions (e.g. a Mott insulator) or Fermi points. These differences completely alter the phase diagram of our model compared to Ref. [25].
The phase diagram for various values of N is summarized in Fig.1. For N > 1, we find that when the kinetic energy of the aforementioned Ising spins is small compared to the kinetic energy of the fermions, the ground state resembles deconfined phase of the Z 2 gauge theory coupled to 2N Dirac fermions (the N = 1 case is special in that the deconfined fermions are unstable to infinitesimal interactions). However, as we will discuss in detail, there is an important distinction between this phase and the deconfined phase of a conventional Z 2 gauge matter theory (Ref. [27][28][29]). Namely, the local Z 2 invariance in our model is an actual symmetry of the Hamiltonian in contrast to the so-called 'gauge symmetry' in gaugematter theories which just corresponds to redundancy in the description of the physical states. This is because we do not project the Hilbert space of our Hamiltonian to gauge invariant states, that is, we do not impose the Gauss's law. Instead, in our model, the Gauss's law in the ground state emerges due to spontaneous symmetry breaking as the locally conserved operators order in a certain definite pattern. One physical significance of this is that in contrast to regular Z 2 gauge theory, in our model, equal space, unequal time non-gauge invariant correlation functions do not vanish identically. Further, there is a finite temperature Ising transition corresponding to the restoration of the aforementioned symmetry.
As the kinetic energy of the Ising degrees of freedom is increased, we find that the aforementioned gapless Dirac deconfined phase enters a symmetry broken confined phase. The nature of symmetry breaking depends on the value of N. For N = 1 it corresponds to a charge density wave, while for N = 2, due to an enlarged symmetry, the symmetry broken phase can correspond to a superconductor or a spin density wave, depending on the sign of infinitesimal symmetry breaking field. The N = 3 case is even richer. We find that after exiting the deconfined phase, the model enters a valence bond solid (VBS) phase, and as the tuning parameter controlling the kinetic energy of Ising spins is tuned further, the valence-bond solid quantum melts and yields a antiferromagnet phase. The phase transition between the VBS and Néel phase, if second-order, will correspond to a deconfined quantum critical point [30], reported earlier in the context of quantum magnets [31,32]. Our results are consistent with a second order transition.

II. MODEL AND SYMMETRIES
Our model reads: Here,ĉ † i i i,α creates a fermion of flavor α in a Wannier state centered around lattice site i i i of a square lattice and the sum runs over bonds of nearest neighbors. Each bond, b b b = i i i, j j j , accommodates a two level system spanned by the Hilbert space |0 b b b and |1 b b b and in this basis,Ẑ i i i, j j j , andX i i i, j j j , corresponds to the Pauli spin matrices 1 0 0 −1 and 0 1 1 0 respectively. h tunes the strength of the transverse field. Note that we have set the coefficient of the hopping term to unity so that all energy scales, including the size of h and the temperature T , will be measured relative to the hopping term.

A. Discrete Symmetries
The model (Eq.1) leads to a number of symmetries. The particle-hole symmetry we consider iŝ and it satisfies P α ,Ĥ = 0.
As an explicit operator,P α = ∏ i (c i i i,α e iK K K·i i i/2 +ĉ † i i i,α e −iK K K·i i i/2 ) where K K K = (π, π). Thus it is a fermionic (bosonic) operator if the total number of lattice sites is odd (even). As it will become clearer later in the paper, it is helpful to define a Z 2 valued order-parameterp i i i that captures the breaking of the particle-hole symmetry, This is nothing but the parity of fermion particle number at site i i i. Under particle-hole transformation,p i i i → −p i i i . A more interesting symmetry corresponds to the local conservation laws: Despite an extensive number of conservation laws, our Hamiltonian is not integrable. This is because the total number of degrees of freedom, N do f , after taking into account L x L y number of conservedQ i i i 's on a L x × L y lattice, are still extensive: N do f = Number of Ising spins Crucially, the operatorsQ i i i anti-commute with the generator of particle-hole symmetry:P −1 αQi i iPα = −Q i i i . Therefore, the low energy theory of our model cannot contain terms that are proportional toQ i i i unless the particle-hole symmetry is spontaneously broken. The anti-commutation also implies that the energy eigenstates cannot simultaneously be particle-hole symmetric and eigenstates ofQ i i i . In particular, particle-hole symmetric eigenstates satisfy Q i i i = 0. However, it's worth emphasizing that even if an energy eigenstate is not particlehole symmetric, it does not necessarily break the particle-hole symmetry spontaneously. For example, one can in principle construct a state |ψ which is an eigenstate of allQ i i i 's, (and therefore cannot be eigenstate ofP α for any α), but nevertheless satisfies ψ|p i i i |ψ = 0, and ψ|p i i ip j j j |ψ → 0 as Another consequence of anticommutation betweenQ i i i and P α is that the spectra of the Hamiltonian is at least doubly degenerate. In fact, on a lattice with an odd number of total sites, the Hamiltonian possesses a N = 2 SUSY [33]. This is because, as discussed inP α is a fermionic operator and together with bosonic operatorQ i i i can be used to define a fermionic SUSY generator Q = H/2P α (1 +Q i i i ) for any α and i i i, so that {Q, Q † } = 2H, Q 2 = 0, and [Q, H] = 0.
At this point it is important to pause and note that even though our Hamiltonian (Eqn.1) bears a strong resemblance to a Z 2 gauge theory coupled to matter fields [27][28][29], ours is actually not a gauge theory due to an important difference. Unlike a matter-gauge theory, we do not impose the Gauss's law constraint at any site. In a conventional Z 2 gauge theory, the fieldX i i i,i i i+n would correspond to the Z 2 electric field along directionn at site i i i, and thus the Gauss's law ('∇ ∇ ∇.E E E = ρ') reads:X i i i,i i i+a a a xX i i i,i i i−a a a xX i i i,i i i+a a a yX i i i,i i i−a a a y (≡ ∇ ∇ ∇.E E E) =p i i i (≡ ρ). Satisfying this operator equation would thus require projecting the Hamiltonian to a specific set of eigenvalues of operatorsQ i i i (namely,Q i i i = 1 [34]). In contrast, we do not impose any constraint onQ i i i 's, which therefore behave like classical degrees of freedom, whose values are allowed to vary from one eigenstate to another, and for a given many-body eigenstate these values will be determined by the intrinsic dynamics of our Hamiltonian.
The local conservation ofQ i i i has important consequences.
or equivalently Thus, a bare single fermionĉ i i i,α as well as the Ising opera-torẐ, are localized along the real space axis, but can propagate along the imaginary time direction. Propagation along the imaginary time implies transitions between differentQ i i i sectors which is not prohibited in our model. As already mentioned, if one were to instead impose a constraint on the values ofQ i i i , it would prohibit propagation along the time direction as well, leading to a local Z 2 gauge redundancy. On the other hand, quantities such as where one attaches a string ofẐ-operators along certain chosen bonds connecting i i i 1 to i i i n , as well as closed paths ofẐoperators, or the expectation value ofX i i i, j j j are all quantities which do not vanish due to symmetry arguments. Again, since the Hilbert space is not restricted to a single choice ofQ i the strings are continuous in space but can be discontinuous in the imaginary time.

B. Continuous Symmetries
The model is clearly invariant under global U(N) rotations in flavor space:ĉ † where U corresponds to an U(N) matrix. We recall that U(N) where the SU(N) component corresponds to the rotations among the N-flavors and the U(1) corresponds to the global particle number conservation c i i i,α → e iθ c i i i,α .
In fact, the full set of continuous symmetries is bigger than U(N) and corresponds to O(2N). This symmetry is manifest if one makes a unitary transformation c j j j,α → ic j j j,α on only one of the sublattices of the square lattice, and then writes the complex fermions in terms of their Majorana components: c j j j,α = γ j j j,α,1 + iγ j j j,α,2 /2. The Hamiltonian then becomes: Since we are in two spatial dimensions, the aforementioned global continuous symmetry can break spontaneously at zero temperature. In many of the examples we will discuss, the symmetry breaking will be manifest in long-ranged correlation functions such as: are the generators of the subgroup SU(N) of the aforementioned O(2N) symmetry. Before moving on, we note that the discrete particle-hole symmetry provides exact mapping between seemingly different correlations functions. For example, let us consider the transverse (α = β ) spin operator Hence, the particle-hole condensate at q q q = K K K = (π, π) is equivalent to a particle-particle condensate at q q q = (0, 0). Equivalently, a charge density wave at q q q = (π, π) transforms to a diagonal spin-density wave at the same wave vector.

III. EFFECTIVE MODEL AT WEAK AND STRONG COUPLING
A. Strong transverse field limit In this limit, we can understand the model in terms of a Su-Schrieffer-Heeger coupling [35] to an Einstein phonon mode of frequency h, but with a truncated phonon Hilbert space. This reading of the model is certainly valid in the high frequency limit where phonon modes are not highly occupied. In particular, and in the |± = 1 √ 2 (|1 + |0 ) basis, we can represent the Ising spins in terms of hard core bosons: At high frequencies, the occupation of the bosonic level will be small so that one can safely relax the hard-core boson con- = 0, and promote the bosonic operators to soft core ones. In this approximation, the analogy to phonons is exact. Integrating out the phonons leads to a retarded interaction resulting in the action: the local hopping term, and the bosonic propagator for −β < τ < β . In the high frequency limit and for constant values of β the interaction becomes instantaneous lim h→∞ D(τ) = δ (τ) and the full Hamiltonian reduces toĤ This interaction was recently considered in the context of SU(N) fermions on the honeycomb lattice in Ref. [36]. The form in Eq. 19 is easy to understand. In the strong coupling limit, the Ising spins are polarized along the x-direction, and fermion hopping is prohibited since this implies flipping of the spin against the transverse field. Virtual processes can nevertheless occur and generate a super exchange energy scale J ∝ 1 2Nh . It is instructive to consider the nature of phases in this limit for a few values of N. First, at large N, the HamiltonianĤ ∞ in Eq.19 can be solved exactly (Ref. [37]) and the ground state corresponds to a valence bond solid state. For N = 1 case, the Hamiltonian is simply H = 1 2Nh ∑ i i i, j j j n i i i − 1 2 n j j j − 1 2 and thus the ground state corresponds to a charge-density wave. For N = 2 case, one finds HereŜ S S i i i is the spin-1/2 operator andη η η i i i =P −1 ↑Ŝ S S i i iP↑ is the 'Anderson's pseudospin' [38], [39]. In the strong coupling limit charge is localized and the Ising spins are frozen in the xdirection wave functions. Such functions have well defined quantum numbersQ i i i . If one chosesQ i i i = −1, then the ground state will be a Néel state for the fermions while the Ising spins point along the x-direction. It is interesting to note that at h = ∞ the kinetic energy vanishes and that all real space charge configurations and thereby allQ i i i sectors are degenerate.

B. Weak transverse field limit
In the h → 0 limit, the imaginary time scale on which the Ising spins can fluctuate diverges such that the Ising spin degrees of freedom become classical. The model then reduces toĤ and one has to find an arrangement of Ising spins which will minimize the energy. Local Z 2 invariance implies that the energy of an eigenstate depends only on the Z 2 fluxF i i i through a plaquette:F Lieb's work [40] implies that the minimal energy is realized by π-flux configurations i.e.F i i i = −1. This means that the coupling to the fermions effectively generates a Z 2 magnetic field term in the Hamiltonian. On the torus with L x × L y sites, there are 2 L x L y configurations of the Ising field that correspond to the πflux through every plaquette, and therefore, an equal number of ground states at h = 0. This huge ground state degeneracy can also be thought of as a result of the conserved quantitieŝ Q i i i that relate different ground states. To see the effect of πflux on fermions, let us chose a two orbital unit cell with unit vectors a a a ± = a a a x ±a a a y such that the intra-unit cell Ising field is −1 while all others are set to 1. As mentioned above, the single particle spectral function depends on the specific choice of Ising fields. On the other hand, the two particle Green functions in the particle-hole channel are independent of the background Ising configuration. In Fig. 2, we plot the ground state dynamical spin-spin correlations function, S(q q q, ω), valid for all N. One observes a particle-hole continuum with gapless excitations at wave vectors connecting the Dirac cones q q q = (π, 0) and q q q = (π, π). One key interest of the QMC simulations will be to study the fate of this particle-hole continuum at non-zero values of h. As one might expect, and in accord with the third law of thermodynamics, quantum fluctuations will lift the zero temperature finite entropy density. In general, at non-zero values of h terms which do not violate symmetries of the model will be dynamically generated. In particular, alongside the magnetic field term of Eq. 23 one can write down, which is invariant under the particle-hole symmetry that sendŝ Q i i i → −Q i i i . Since this is nothing but a classical Ising Hamiltonian, one can foresee that there will likely be a finite temperature Ising transition below which theQ i i i 's spontaneously order. As discussed in the following sections, this symmetry breaking will result in an emergent Z 2 gauge structure in the ground state phase diagram and would also determine the nature of both confined and deconfined phases.

IV. METHODOLOGY
In the Majorana representation of Eq. 11 the O(2N) symmetry of our model is manifest. As a consequence, sign free quantum Monte Carlo (QMC) simulations can be carried out for all values of N [41]. We have used a standard implementation of the finite temperature auxiliary field approach [15,16,42] where the trace runs over the fermionic and Ising degrees of freedom and β corresponds to the inverse temperature. For Hubbard type model simulations, the auxiliary field corresponds to a Hubbard-Stratonovitch field with no explicit imaginary time dynamics [15,16,42]. Here, in contrast, the auxiliary field corresponds to the z-component of the Ising field, and its imaginary time dynamics is controlled by the transverse field. We have adopted a single spin flip algorithm, which turns out to be efficient (high acceptance rates) at large values of h where the field oscillates rapidly in imaginary time. At weak couplings, the Ising field becomes classical and the single spin-flip update in space-time becomes increasingly inefficient. This renders simulations at low values of h expensive. Unless mentioned otherwise we have used a Trotter time step ∆τ = 0.2 so as to allow for a bond decomposition of the infinitesimal time propagator: corresponds to the Hamiltonian on a single bond. Within the auxiliary field QMC method one can compute equal time as well as imaginary time displaced correlation functions. We have carried out the analytical continuation with a stochastic implementation of the Maximum Entropy method [43,44]. Eq. 7, which states that the single-particle Green's function is independent of momentum, provides an ideal test for the QMC as well as for the analytical continuation since the locality of the single particle Green function results from the Monte Carlo sampling. Fig. 3 shows the single particle spectral function A(k k k, ω) = −Im G ret (k k k, ω) for N = 1 model which is obtained from Wick rotation of the imaginary time Green function G(k k k, τ) = − Tĉ k k k (τ)ĉ † k k k (0) . The k k k-independence of the data should be seen as a measure of our spectral resolution.
FIG. 3. (color online) Single-particle spectral function A(k k k, ω) for the N = 1 model. As argued in the text, A(k k k, ω) is k k k-independent.

V. NUMERICAL RESULTS
In this section we will describe our numerical results obtained from finite temperature auxiliary field simulations.

Summary of Results
A sketch of the zero temperature phase diagram in the transverse field h versus N plane is shown in Fig. 1. We will argue below that for all values of N, there is a finite temperature Ising transition at which theQ i i i operators order. Since Q i i i ,Ĥ = 0, this is akin to spontaneous symmetry breaking in a classical Ising model. For even (odd) values of N we observe (anti) ferromagnetic ordering of theQ i i i 's so that at zero temperature we can understand our results in terms of an effective Z 2 gauge theory coupled to 2N Dirac fermions, with the understanding that non-gauge invariant unequal time, equalspace correlators do not vanish, as discussed earlier. The (anti) ferromagnetic ordering for even (odd) values of N is consistent with the fact that perturbatively, the couplings K i, j in Eq. 24 for nearest neighbour vertices are proportional to (−1) N+1 .
Spontaneous breaking ofQ i i i → −Q i i i generates term proportional to ∑ i e iNK K K·i i ip i i i in the effective Hamiltonian, leading to spontaneous charge ordering of fermions below the finite temperature transition. The relevance/irrelevance of this term determines the fate of h = 0 gapless Dirac fermions (Fig. 2) at infinitesimal h. When N = 1, this term is a fermion bilinear proportional to ∑ i e iK K K·i i i (1 − 2n i i i ) and thus leads to spontaneous charge-ordering and mass gap for fermions at infinitesimal h. For N = 2, this term is proportional to the Hubbard term where the overall sign of this term is determined by whether Q i i i = 1 or −1. The Hubbard term is irrelevant for gapless Dirac fermions, and therefore one expects to find a stable Z 2 Dirac deconfined phase (Z 2 D in Fig.1). Similarly, when N = 3, the correspond- 3 ) is irrelevant at the h = 0 point, leading again to a stable Z 2 D phase.
As h is increased, the aforementioned Dirac deconfined phase undergoes a transition to a symmetry broken confined phase whose nature depends on the value of N (Fig.1). For N = 2, we find a Néel antiferromagnet, or a superconductor depending on whether Q i i i = −1 or 1, respectively. For N = 3, we find two symmetry broken phases, a dimerized (≡ Valence Bond Solid) phase separated from the Z 2 D phase again by a seemingly continuous transition, and a Néel phase at larger values of h. The transition between the spin-dimerized state and the Néel phase is apparently also continuous and thus provides an example of deconfined quantum critical point [30].

A. N=1
As already discussed briefly in Sec.II, in the large-h limit, our model Hamiltonian (Eq.1) maps onto (Eq. 19) so that the ground state is a charge density wave. In this limit X b b b 1 and the charge ordering follows a K K K = (π, π) modulation. Consequently theQ i i i 's, by definition (Eq.6), order antiferromagnetically. Therefore, as a function of temperature, there must be a finite temperature transition in the 2D Ising universality class below which theQ i i i develop long range order. In the absence of a particle-hole symmetry breaking field, Q i i i = 0 on any finite lattice and one has to measure correlation functions to detect ordering. The Q i i iQ j j j correlation functions turn out to be very noisy and we have found it more convenient to include a small symmetry breaking field: Such a term in Hamiltonian does not introduce a negative sign problem [45] and pins the CDW. Since Q i i i ,Ĥ s = 0,Q i i i remains a good quantum number at finite value of h s . In fact H s acts on the Ising model of theQ i i i variables as a staggered magnetic field. Figure 4(a) plots at a small value of the particle-hole symmetry breaking field, h s = 0.01. At low temperatures we find that Q K K K scales to unity for all values of h. SinceQ i i i is a good quantum number, even in the presence of the symmetry breaking field, this means that the ground state haŝ and the gauge constraint is dynamically generated. The temperature scale at which Q K K K takes half of its maximal value defines T 1/2   The product of these two terms, which will also be generated, is proportional to ∑ i i i e iK K K·i i ip i i i and leads to charge-ordering. To confirm this, we have hence computed As apparent from Fig. 4(b) the energy scale at which p K K K picks up, matches T 1/2 defined in Eq.30. As a further test we have computed charge-charge correlation functions in the absence of symmetry breaking field: N(q q q) = 1 4 ∑ r r r e iq q q·r r r p r r rp0 0 0 (32) At low temperatures, ordering occurs at the antiferromagnetic wave vector K K K = (π, π) and as shown in Fig. 4(c) a reasonable data collapse is obtained using the 2D Ising exponents, η = 1/4 and ν = 1. For the value of h = 1 considered in Fig. 4(d) we obtain T c 0.17 which stands, to a first approximation, in agreement with the value of T 1/2 listed in Fig. 4(b). Fig. 5 plots 'gauge invariant' quantities X b b b and F i i i (Eq.22) corresponding to the Ising field. We note that the av-  erage electric field is given by where F is the free energy, and therefore provides a direct access to the behavior of free energy as a function of h. At small values of h one observes strong finite size and temperature effects for both quantities. Our lowest temperature results support a smooth crossover from This implies that there is no level crossing associated with a change of ground state quantum numbers and that the gauge constraint remains enforced for all values of the transverse field. At h = 0 all plaquettes are threaded by a π-flux static magnetic field and no electric field is present. Starting form this limit gauge fluctuations introduce pairs of Z 2 vortices -which are nothing but plaquettes with zero flux -such that as a function of h the magnetic field grows. The fluctuations of the magnetic field induce the onset of the electric field.
We now discuss zero-temperature phase diagram. The CDW order parameter, at low temperatures and as a function of h is plotted in Fig. 6. Consistent with the fact that at h = 0, the spectrum is gapless (Fig. 2), and correspondingly T 1/2 → 0 as h → 0 (Fig. 4), the CDW order parameter also approaches zero as h → 0. We expect that Q i i i = 0 for any non-zero h, and therefore the Z 2 D phase is stable only at h = 0.
In Fig. 7 we consider the dynamical properties of the model at h = 0.1 and β = 40. This choice of parameters places us below T 1/2 and we will discuss the data in terms of the ground state with quantum numberŝ  The dynamical charge structure factor reads, withN The intermediate states, |Ψ n , have the sameQ i i i quantum numbers as in the ground state since n i i i ,Q i i i = 0. Thereby, the dynamical charge structure factor of Fig. 7 is identical to that one would obtain when imposing exactly the gauge constraint on the Hilbert space. Within this constrained Hilbert space, the elementary excitations are closed Wilson loops of Z b b b operators in Euclidian space-time as well as open strings ofẐ b b b operators with particle-hole excitations at the extremities again in 2 + 1 dimensions. For the considered values of h the ground state shows long range CDW order and the gauge string binds particle-hole excitations. This corresponds to the confined phase of the gauge field and leads to the central peak feature at the wave vector K K K = (π, π) apparent in Fig. 7. In the limit h → ∞ this becomes the dominant feature of the spectral function. At higher energies or on shorter time scales, the string does not bind particle-hole excitations and a continuum very reminiscent of the h = 0 limit is apparent in the spectral function.

B. N=2
As discussed in the Sec.II B, for N = 2 the global symmetry is O(4). To understand the physical content of this symmetry, it is instructive to consider the limit h → ∞: Here,Ŝ S S i i i = 1 2 ∑ σ ,σ ĉ † i i i,σ σ σ σ σ ,σ ĉ i i i,σ corresponds to a fermionic representation of the spin operator andη η η i i i =P −1 1Ŝ S S i i iP1 whereP 1 corresponds to particle-hole symmetry transformation defined in Eq. 2. Since the η η η-operators are derived from the spin operators with a similarity transformation they satisfy the SU(2) algebra as well. Furthermore, η j j j,α ,Ŝ i i i,β = 0 so that at the level of Lie algebra, the O(4) symmetry can be understood as O(4) ∼ = SU S (2) ⊗ SU η (2) [46], [47].
The Ising order parameter corresponding to particle-hole symmetry breaking (Eq. 4) relates to theŜ S Sandη η η-operators viap If particle-hole symmetry breaks, then there are two possibilities. If p i i i < 0, then the SU S (2) ⊗ SU η (2) is reduced to SU S (2). In the extreme limit of p = −1, sites are singly occupied, and the model reduces to a spin-only model. Away from this limit, there are charge fluctuations per site, but the universal properties (such as the phase transition out of the deconfined phase), which depends only on symmetries, can still be understood in terms of only spin fluctuations. The other possibility is that p i i i > 0, in which case the sites are predominantly empty or doubly occupied, leading to a charge-only model in the extreme limit p = 1. In this case, the underlying symmetry is SU η (2) and the phase diagram can be understood in terms of fluctuations of a three component vector consisting of CDW (one component) and s-wave superconductor ( two components, since it has both real and imaginary components). Since the superconducting phase, whenever it exists, is always degenerate with CDW due to SU η (2) symmetry, henceforth we will call the corresponding phase as just superconductor.
To detect the possibility of Ising transitions corresponding to ordering ofQ i i i andp i i i we compute the temperature dependence of the internal energy for a range of values of h (Fig. 8).   Note that in contrast to theQ i i i 's,p i i i is not a good quantum number since p i i i ,Ĥ = 0. The Ising transition in two-dimensions has a specific heat exponent α = 0 which implies that the specific heat is logarithmic divergent at the transition temperature. This results in an infinite slope of the internal energy at T c and in the thermodynamic limit. An arrow in each plot points to the temperature where singularities are apparent. It is beyond the scope of our calculations to pin down the details of the singularity. In particular, when T c is small very large systems sizes are required to resolve the two dimensional nature of the singularity. Within our resolution we observe only one transition which we will now argue corresponds to the ordering of theQ i i i 's as well asp i i i , that is, breaking of particle-hole symmetry.
To explicitly see the ordering ofQ i i i andp i i i , we add a small attractive Hubbard term to the model,Ĥ U = −U ∑ i i i (n i i i,1 − 1/2) (n i i i,2 − 1/2). This does not introduce a sign problem [48], but breaks particle-hole symmetry and effectively acts as a uniform magnetic field for theQ i i i variables, alongside the Ising interaction term of Eq. 24. Note that Ĥ U ,Q i i i = 0 so thatQ i i i remains a good quantum number upon inclusion of this symmetry breaking term. Fig. 9 shows at low temperature implies ferromagnetic coupling between theQ i i i variables in Eq. 24. Since the Hubbard term acts as a Zeeman field for particle-hole symmetry as well, we can compute p i i i to check if particle-hole symmetry is broken spontaneously at the same energy scale. If so, one expects p i i i to show a marked growth at the energy scale where Q i i i ≈ 0.5. Fig. 9(b) supports this. Fig. 9(c) summarizes our estimate of T c as obtained from the internal energy and by explicitly including a symmetry breaking term.
Below T c , any change in theQ i i i quantum numbers as a function of h should correspond to a level crossing. To rule out this scenario, we plot the average electric (see Eq. 33) and magnetic field per plaquette (Eq. 22) as a function of h. Size and temperature effects are strong in the weak coupling limit where T c is small. However, our lowest temperature results on the L = 16 lattice are consistent with a smooth evolution of the magnetic flux and electric field. Hence, at an energy scale set by T c the constraint is dynamically generated such that at T = 0 we can interpret our findings in terms of a Z 2 lattice gauge theory withQ i i i = 1, orQ i i i = −1 coupled to fermions.
As discussed above, the pattern of ordering forQ i i i determines the pattern of symmetry breaking of the continuous global symmetries in the confined phase. Specifically, superconducting phase corresponds toQ i i i = 1 and p i i i > 0 while the SDW phase corresponds toQ i i i = −1 and p i i i < 0. One can in fact explicitly map a superconducting state to an SDW state by particle-hole transformation. Consider a superconducting ground state |Ψ , so that Since P ↑ ,Ĥ U = 0 and P ↑ ,Q i i i = 0, it follows that one can generate an orthogonal state |Ψ =P ↑ |Ψ which satisfies If the state |Ψ corresponds to a superconductor, then |Ψ corresponds to SDW, since Ψ|Ŝ S S i i i .Ŝ S S j j j |Ψ = Ψ |η η η i i i .η η η j j j |Ψ for all i i i, j j j. In the limit U → 0, |Ψ would also correspond to an eigenstate of the Hamiltonian. As an aside, as discussed in the introduction, when the number of sites is odd, then |Ψ and |Ψ are supersymmetric partners [33], and thus N = 2 supersymmetry relates superconducting and SDW order parameters. Similar dichotomy between superconductor and SDW was also explored in [49]. Let us investigate the zero temperature in the language of spins i.e. in the sectorQ i i i = −1. We compute the spin-spin correlations given by: The local moment of the antiferromagnetic order is then defined as: Note that Ŝ α β (i i i),Q i i i = 0 such that spin excitations do not violate the dynamically generated constraint at T = 0. Fig. 11 shows the finite size scaling of m for various value of h. We see that extrapolation to the thermodynamic limit with a second order form in 1/L yields a magnetization which vanishes at h 0.1. Below h 0.1 we have tested for bond dimerized orders but found no positive signal. In this regime, we will interpret our findings in terms of a Z 2 -Dirac deconfined phase of our effective lattice gauge theory. At T = 0 whereQ i i i = −1, the gauge invariant elementary excitations are pairs of spinons tied by a gauge field string ofẐ b b b operators as well as closed loops ofẐ b b b 's. Below h 0.1 the gauge field is deconfined and the magnetic flux per plaquette is close to π. In particular one expects the dynamical spin structure factor to bear similarities with the particle-hole continuum of Dirac fermions shown in Fig. 2. In fact the Fig. 2 would correspond to the mean-field theory of the Z 2 -Dirac deconfined phase. Fig. 12 plots the dynamical spin-spin correlation functions as a function of h. As apparent, in the small h limit we indeed see a continuum of excitations with soft features at wave vectors (0, π) and (π, π) corresponding to transitions between Dirac nodes. As h grows spectral weight accumulates around K K K = (π, π) and a distinct spin-density wave feature emerges. As in the N = 1 case, in the proximity of the transition the continuum of fractionalized excitations is visible at high energies whereas the low energy sector is dominated by the Goldstone mode.
To see the difference between the small h vs large h phase, we also considered the spin-susceptibility, Some care has to be taken when considering this quantity, since Monte Carlo averages over differentQ i i i sectors, and it is only in the low temperature limit that this quantity will correspond to that of the lattice gauge theory. Our results are shown in Fig. 13. As apparent, the two phases have very distinct low energy signatures. Taking into account finite size effects, the Z 2 -Dirac deconfined phase is characterized by a χ s ∝ T in accord with the mean-field theory of this state. On the other hand, the confined phase is consistent with a constant low temperature spin-susceptibility originating from spin-wave excitations.

C. N=3
We will next discuss the N = 3 model (Eq.1) where alongside a deconfined phase at small h and magnetically ordered phase at large h, a spin dimerized state emerges at intermediate values of h. As in the N = 1 case, theQ i i i 's order antiferromagnetically such that we can pick up the ordering temperature by implementing a small staggered field as in Eq. 27. Our results are plotted in Fig. 14(a) and (b). As expected, anti-ferromagneitc ordering ofQ i i i andp i i i occurs at the same energy scale. Being a constant of motion even in the presence of the small staggered chemical potentialQ K K K (see Eq. 28) saturates to unity. As a technical aside, we note for values of h < 1/3 the Monte Carlo autocorrelation time for aforementioned quantities becomes large thus rendering the calculation difficult.
In Fig. 14(c) and (d) we plot the magnetic flux as well as the electric field as a function of h. Size and temperature ef- fects are large at a small values of h. Nevertheless, our data is consistent with absence of level crossing in the ground state quantum numbers since owing to Eq. 33 this would result in a jump in the electric field X b b b . The anti-ferromagnetic ordering ofp i i i along with the constraint that there are on average one and a half electrons per unit cell, gives us two options for the resulting charge ordering: I. One can place two charges on sub-lattice A and one charge on sub-lattice B. This choice leaves the SU(3) color degree of freedom unquenched. In particular, on sub-lattice A, two SU(3) spins combine to give 3 ⊗ 3 = 3 ⊕ 6 where 3 corresponds to the fundamental representation, 3 to the antisymmetric complex conjugate representation with Young tableau  [50]. For lower values of h we again use a polynomial fit and impose a zero intercept with the y-axis. In Fig. (b) we plot the extrapolated value of the magnetization as well as of the VBS order parameter. At h = 3/8 both VBS and SDW order parameters are consistent with a power-law decay.
consisting of a single column and two rows and 6 to the symmetric representation with Young tableau consisting of a single row and two columns. Since we are working with fermions, only the 3 representation is relevant and has states given by Here summation over repeated indices, running over the three colors of SU(3), is implied. Neglecting the site index i i i and assuming that under an infinitesimal SU(3) rotationĉ † γ → (1 + iεe e e ·T T T ) γ,δĉ † δ with T T T , the vector of Gell-Mann matrices, e e e a unit vector in R 8 , and ε a real infinitesimal, one will show that ∆ † α → 1 − iεe e e ·T T T α,δ ∆ † δ . On sub-lattice B the single electron corresponds to the fundamental representation of SU (3). With this assignment of representations, Néel states, as well spin-dimerized states can naturally occur. For the Néel state we break the SU(3) symmetry by arbitrarily choosing a color α so as to define the state: Fluctuations around this state will produce spin waves. On the other hand, a possible VBS state with q q q = (π, 0) order is given by: This state remains invariant under SU(3) rotations and thereby defines a singlet. II. The other charge ordering pattern consistent with the antiferromagnetic structure of thep i i i 's, corresponds to three particles on sub-lattice B and none on sub-lattice A. This assignment quenches the spin degree of freedom. One might expect that SDW states which possess low lying Goldstone modes are energetically favorable. Our QMC results are consistent with this expectation. In Fig.15 we plot the equal time spin-spin correlation and corresponding staggered moment m. Antiferromagnetic spin ordering is present down to h 3/8. [51] At lower values of h, the data is consistent with the vanishing of SDW.
The uniform spin susceptibility plotted in Fig. 16 gives a very clear account on the nature of the phases one expects as a function of h. At h = 1 and at low temperatures, the data is consistent with the SDW state with constant density of states originating from the Goldstone modes. At h = 1/12 we observe, similar to N = 2 case at h = 0.075, the characteristic spin susceptibility of a Dirac deconfined phase. There is however a distinct intermediate phase possessing no low lying spin excitations. As hinted above, a very natural candidate to account for this is spin-dimerization which breaks the C 4 lattice symmetry. To verify this, we compute the dimer-dimer correlation functions: [52]. From the above correlation function, we define the dimer order parameter as: Fig. 17 shows a finite size scaling of the above VBS order parameter at wave vector q q q = (π, 0). Extrapolating to the thermodynamic limit necessarily implies a choice for the fit function. In the parameter range where we observe no low lying spin-excitations in the susceptibility, we use a polynomial fit up to second order in 1/L. This fit indeed confirms long range valence-bond order in the vicinity of h = 1/6.
In Fig. 18 we show the temperature dependence of the VBS and SDW correlation functions at q q q = (π, 0) and at q q q = (π, π) respectively and in the Z 2 D, VBS and SDW phases. The VBS order breaks the lattice C 4 symmetry and therefore one correspondingly expects a finite temperature Kosterlitz-Thouless transition below which S VBS (π, 0) diverges as a function of system size. The SDW ordering breaks a continuous symmetry and at finite temperatures the real space spin-spin correlations will decay exponentially if the total system size is bigger than the correlation length. In our model, there is also the Ising transition of the Q Q Q i i i 's which will feed into the spin and dimer correlation functions. At h = 1 we can evaluate T 1/2 , Q K K K (T 1/2 ) = 1/2, from the data shown in Fig. 14(a). This energy scales, which tracks the Ising transition temperature, is shown in Fig. 18 (c). Here, the VBS correlations show strong size effects above the T 1/2 scale suggesting that the Kosterlitz-Thouless transition temperature exceeds the ordering temperature of the Q Q Q i i i 's. At T 1/2 charge ordering equally occurs (see Fig. 14(b)) leading to a dynamical selection of the representations of SU(3) on each sub-lattice. The fact that at T 1/2 the antiferromagnetic spin correlations (Fig. 18 (f)) grow agrees with the 3 and3 representation pattern. The drop in the VBS order parameter at this energy scale illustrates the competition between these two ordered states and the data at h = 1 supports dominant SDW correlations in the low temperature limit. We note that the data of Fig. 18 (f) seemingly supports long range order at finite temperature. However, recall that the antiferromagnetic correlation length grows exponentially with inverse temperature [53] so that when it exceeds our largest lattice size the QMC data falsely suggests long range ordering. Our results at h = 1/6 in the VBS phase are plotted in Fig. 18 (b),(e). Comparing with the h=1 data set, one observes very similar finite temperature behavior and we again interpret the sharp downturn in the VBS correlation functions as a measure of the Ising transition. However, and in contrast to the h = 1 data set, the low energy results are consistent with long range dimer order. In the Z 2 D phase (see Fig 18 (a),(d)) the Ising transition is again manifest in the sharp upturn in both the dimer and the spin correlations. For this value of the coupling both spin and dimer correlation grow below the Ising transition.
Overall the data is consistent with a deconfined phase at small values of h, a valence bond solid at intermediate values of h and finally antiferromagnetic order at large values h. The evolution of the order parameters is show in Fig. 15(b). From this we observe that the SU(3) model not only captures a seemingly continuous transition between the Z 2 Dirac deconfined and the VBS, similar to the N = 2 case, but also a transition from the spin-dimerized state (VBS) to the SDW phase. Since we observe no jump in the electric field and the order parameters also evolve continuously, Eq. 33 leads to the conclusion that all transitions are continuous. In particular, the VBS and SDW phases is then expected to be an instance of a deconfined quantum critical point occurring at h 3/8, similar to those found in quantum magnets [32]. The low energy theory for this transition has been argued to be three complex bosons coupled to a non-compact U(1) gauge field (noncompact CP 2 theory) [30]. Within our limited system sizes we were not able to confirm the literature value of the exponents [50,54]. Note however that even for the largest values of L in Ref. [54] the ν exponent has still not converge. We notice that irrespective of the value of ν, β m < β V BS where β corresponds to the order parameter exponent. This inequality between β m and β V BS is seemingly consistent with Fig. 15(b) where VBS order parameter rises more sharply compared to the SDW as one moves away from the transition on the two sides. Fig. 19 plots the evolution of the spin dynamical structure factor as a function of h. As for the SU(2) case, Q i i i ,Ŝ β α ( j j j) = 0 such that intermediate states in the Lehmann representation of the zero temperature spectral function do not violate the dynamically generated constraint. Since we are running through two phase transition in rather close proximity, we have to disentangle features stemming from the Z 2 Dirac deconfined phase and from the deconfined quantum critical point. The characteristic of the Z 2 D phase at h = 1/12 is the spinon continuum, the high energy features of which are very visible across the transition, and in the VBS phase (h = 1/6) but then fade away when entering the SDW phase (h ≥ 3/8) and are completely absent in the strong coupling limit at h = 3. At this coupling strength, one observes a pronounced Goldstone mode in the vicinity of K K K = (π, π) accounting for the spin order. Starting from this point and turning down the transverse field h results in a loss of weight of the spin wave mode. Upon close inspection, one equally sees that it shifts to higher energies.
Another striking feature is the overall broadness of the spectral function in the vicinity of the antiferromagnetic wave vector and in proximity of the VBS to SDW transition at h 1/3. Qualitatively, this is consistent with a deconfined quantum critical point which is characterized by a large value of the anomalous dimension η characterizing a complete breakdown of well defined quasiparticles at criticality. Of course, since the spin susceptibility at criticality is expected to be of the form: with c the spin velocity, any non-zero η implies absence of quasiparticles. However at a conventional Wilson-Fisher fixed point in 2+1 dimensions, which describes a phase transition between an SDW and a featureless paramagnet in the Heisenberg universality class η ≈ 0.036 [55,56], and correspondingly there are essentially no broad features visible in S(q q q, ω) [57]. For the SU(3) J-Q model which shows a deconfined critical point, η 0.4 [50] . We note that the Z 2 D to SDW transition observed in the SU(2) model is equally characterized by a large value anomalous exponent which within the ε-expansion reads η = 0.8. In the vicinity of the transition (see Fig. 12 a) and b)) broad features at the antiferromagnetic wave vector are again apparent. The data set at h = 3 (see Fig. 19e)) corresponds to the strong coupling limit where the spin wave velocity vanishes as 1 2Nh . Upon comparison between the data sets of Fig. 19d) and e) one notices a marked reduction of the spin wave velocity. Another striking feature of the data, is the lack of a sharp spin-wave mode as seen in the SU(2) case at our strongest coupling (see Fig. 12e)). In fact in the vicinity of the q q q = (π, 0), one detects two distinct features. We interpret this in terms of charge fluctuations around the mean-field spin-density wave state of Eq. 46. Setting α = 1 in this equation and applying HamiltonianĤ ∞ will generate spin-flip excitations on bonds with α = 2, 3. This set of excitations lead to the spin-density wave mode. One will however equally generate charge fluctuations where one site of the bond is left empty and the other accommodates three particles. We interpret the additional features observed in the data in terms of these charge fluctuations.
At h = ∞ spin and charge excitations are degenerate and the spin wave velocity as well as the energy cost of the charge excitation are tied to the same scale 1 2Nh such that we cannot dissociate both features. Furthermore one easily sees that such charge fluctuations do not occur in the SU(2) case, where we observe a sharp spin-wave mode (Fig.12e).

VI. SUMMARY AND FUTURE DIRECTIONS
Perhaps the most attractive feature of our model (Eq.1) is that despite its simplicity and the absence of fermion sign problem, the phase diagram is rather rich (Fig.1) and features several highly entangled phases and quantum phase transitions. In particular, in the Z 2 deconfined phase, topological degrees of freedom, namely, the dynamically generated Z 2 gauge field is coupled to gapless fermions. From an entanglement perspective, the Z 2 Dirac phase will result in sharp universal footprints in their entanglement structure and this provides an attractive opportunity to apply some of the recently developed techniques to calculate entanglement entropies and entanglement spectrum [58][59][60][61][62]. Specifically, the universal part of the entanglement in the Z 2 Dirac phase is expected to follow the relation γ Z 2 D = γ D + log(2) where γ D is the universal value for the same shaped region in a regular Dirac phase [63]. It would be similarly interesting to study the entanglement structure at the transition from Z 2 D to symmetry broken phases for N = 2 and N = 3, as well as at the deconfined critical point between SDW and VBS.
The transitions from Z 2 D phase to the SDW/SC for N = 2 and to VBS for N = 3 are rather exotic since they appear to be continuous while involving both symmetry breaking, and confinement of the gauge field [64]. We don't have an understanding of these transitions and it's a worthwhile direction to pursue.
Although our model is not microscopically motivated by any specific material(s), the resulting phenomenology does share several features with spin-liquid candidate materials [65][66][67][68] which also exhibit broad, non-quasiparticle features in their excitation spectrum. It is also worth noting that there have been recent proposals for realizing dynamical gaugematter theories in cold atomic lattices where instead of imposing the gauge constraint as a penalty term, the Hamiltonians such as ours (Eq.1) appear naturally due to angular momentum conservation [69,70].
Our model features an interesting interplay of several symmetries -an extensive number of conserved operatorsQ i i i that anti-commute with the particle-hole symmetry, as well as continuous O(2N) symmetry that allows one to perform sign free Monte Carlo for all N. The particle-hole symmetry is the reason that our model hosts a finite temperature transition whereQ i i i spontaneously orders leading to an effective 'Gauss's law' for a Z 2 gauge-matter theory, unlike a conventional gauge theory [27,29] where the Gauss's law is imposed explicitly on the Hilbert space. Contrast this also with 2D Toric code [14] where the Gauss's law is imposed energetically by explicitly adding a term that is proportional to ∑ i i iQ i i i = ∑ i i iXi i i,i i i+a a a xX i i i,i i i−a a a xX i i i,i i i+a a a yX i i i,i i i−a a a y . Since this term breaks theQ i i i → −Q i i i symmetry explicitly, there is no finite temperature transition in the 2D Toric code. The fieldsQ i i i in our model are of course classical in the sense that they are conserved and don't have any dynamics. It would be desirable to extend the idea of studying deconfined phases of gauge theories without imposing the gauge constraint to other gauge groups as well, particularly to the case of U(1) gauge theory relevant to several spin-liquid candidates. Imposing gauge constraint numerically is often challenging, and our results indicate that at least in certain gauge theories, is not necessarily required.
On the above note, it would also be interesting to explore Toric code [14] like Hamiltonians with inbuiltQ i i i → −Q i i i symmetry, and study the implications of the melting of the Gauss's law at finite temperature. This could be especially interesting in 4D Toric code where there exist finite temperature quantum memory [71]. Is the quantum memory lost as soon as the thermal expectation value Q i i i vanishes?
In our models, the finite temperature transition corresponding to the ordering ofQ i i i coincided with spontaneous breaking of particle-hole symmetry for all values of N. This is because X i i i,i i i+a a a xX i i i,i i i−a a a xX i i i,i i i+a a a yX i i i,i i i−a a a y = 0. It is interesting to contemplate the possiblity of a phase where this is not the case such thatQ i i i order while the spatial correlation functions ofp i i i are short-ranged. Consider for example the model: whereQ i i i =X i i i,i i i+a a a xX i i i,i i i−a a a xX i i i,i i i+a a a yX i i i,i i i−a a a y . One can implement â Q i i i → −Q i i i symmetry in the above Hamiltonian by demanding thatĤ be invariant under the transformationÛ −1ĤÛ whereÛ = ∏ i i i, j j j Ẑ i i i, j j j ∏ i i iPi i i where the prime on the first product denotes that it is taken over all vertical (horizontol) links either only for even rows (columns) or only for odd rows (columns). This transformation takesQ i i i → −Q i i i ,Q i i i →Q i i i andX i i i, j j j → −X i i i, j j j . Therefore, even though Q i i i = 0 in this model, the order parameter for the charge-ordering, namely, p i i i will be non-zero only if theÛ symmetry is spontaneously broken. A phase whereQ i i i order whilep i i i remain disordered would necessarily involve interesting entanglement between the fermions and the Ising degrees of freedom because one requires that X i i i,i i i+a a a xX i i i,i i i−a a a xX i i i,i i i+a a a yX i i i,i i i−a a a yp i i i = 0 while p i i i = 0.
The N = 3 model has the richest phase diagram among the models we study. In particular, apart from the continuous transition between Z 2 D and VBS, the phase diagram contains two distinct symmetry broken phases, the VBS and the SDW phase. Remarkably, we find that these two phases are separated by a second order transition (Fig. 15b), which is thus expected to belong to a non-compact CP 2 universality class [30]. At this deconfined quantum critical point, we also study the imaginary-time dynamics via the calculation of structure factor S(q, ω) (Fig. 19, and found a rather broad continuum of excitations, consistent with the theoretical expectation that at this transition, the effective description is in terms of fractionalized 'spinons'. To our knowledge, the dynamics at the deconfined quantum criticality has not been studied numerically in the past, and given rather distinct features seen in our work, it would be worthwhile to pursue this direction in quantum spin-models where much larger sizes are accessible via bosonic Quantum Monte Carlo such as Stochastic Series Expansion (SSE) [31].
Some cautionary statements are in order for the N = 3 case. We found it challenging to perform scaling collapse at the both quantum phase transitions, despite no indication of a firstorder transition in the behavior of the order parameter. One reason for this could be that the two zero temperature transitions are rather close by in the phase diagram. This is visible in the behavior of the structure factor as h is tuned (Fig. 19) where a diffuse continuum of excitations exists all the way from h = 0 (Z 2 Dirac phase) to h ∼ 1/3 (VBS to SDW transition). It's also worth mentioning that the system sizes we are working (L ≤ 16), are an order of magnitude smaller than the ones for spin-models which exhibit deconfined quantum criticality.
We notice that when N is even, our model does not suffer from sign problem even in the presence of the a non-zero chemical potential on arbitrary lattices. Thus it would be interesting to extend our study to a system with finite density of fermions.
We would like to conclude with a conjecture only tangentially related to the rest of the paper. Consider the Hamiltonian in Eq. 19 for N = 2 on an arbitrary graph,Ĥ ∞ =Ĥ η η η + H σ σ σ whereĤ η η η = ∑ i i i, j j j J i i i, j j j η z i i iη . Assuming J i i i, j j j ≥ 0, the HamiltonianĤ ∞ is sign problem free on an arbitrary graph because one can use Hubbard-Stratanovich decomposition to write the partition function as a sum of nonnegative fermion determinants. On the other hand, considering its decomposition into spin-HamiltoniansĤ η η η andĤ s s s , onlyĤ η η η is sign problem free on an arbitrary graph due to all non-positive off-diagonal elements in a local basis. Since [η η η, S S S] = 0, the simulatability ofĤ ∞ suggests that its ground state is likely identical to that ofĤ η η η . That is, in the ground state one expects no singly occupied sites -half of the sites are unoccupied while the other half are doubly occupied, so that onlyĤ η η η acts non-trivially. This leads to the following conjecture: Given a set of J i j ≥ 0, the ground state energy ofĤ − = ∑ i i i, j j j J i i i, j j j σ z i i i σ z j j j − 1 2 σ + i i i σ − j j j + σ − i i i σ + j j j is less than or equal to the ground state energy ofĤ + = ∑ i i i, j j j J i i i, j j j σ z i i i σ z j j j + 1 2 σ + i i i σ − j j j + σ − i i i σ + j j j . Clearly, on a bipartite graph this conjectured inequality is always saturated. To test this conjecture, we performed exact diagonalization study for systems upto 10 qubits when J i j are chosen randomly from a uniform distribution on [0, 1] for a sample size of 1000. We didn't find any counterexample.
Prior to submission of this paper, we were aware of related work by Snir Gazit, Mohit Randeria and Ashvin Vishwanath [72]. We thank them for sharing their draft.