Fermi-surface reconstruction without symmetry breaking
Snir Gazit Fakher F. Assaad Subir Sachdev


	
	
	Racah Institute of Physics and The Fritz Haber Research Center for Molecular Dynamics, The Hebrew University of Jerusalem, Jerusalem 91904, Israel
	
	
	Institut für Theoretische Physik und Astrophysik,
		Universität Würzburg, Am Hubland, D-97074 Würzburg, Germany
	Würzburg-Dresden Cluster of Excellence ct.qmat,
		Universität Würzburg, Am Hubland, D-97074 Würzburg, Germany
	
	
	Department of Physics, Harvard University, Cambridge MA 02138, USA
	
		We present a sign-problem free quantum Monte Carlo study of a model that exhibits quantum phase transitions without symmetry breaking and associated changes in the size of the Fermi surface. The model is an Ising gauge theory on the square lattice coupled to an Ising matter field and spinful `orthogonal' fermions at half-filling, both carrying Ising gauge charges. In contrast to previous studies, our model hosts an electron-like, gauge-neutral fermion excitation providing access to Fermi liquid phases. One of the phases of the model is a previously studied orthogonal semi-metal, which has ℤ_2 topological order, and Luttinger-volume violating Fermi points with gapless orthogonal fermion excitations. We elucidate the global phase diagram of the model: along with a conventional Fermi liquid phase with a large Luttinger-volume Fermi surface, we also find a 
			`deconfined' Fermi liquid in which the large Fermi surface co-exists with fractionalized excitations. We present results for the electron spectral function, showing its evolution from the orthogonal semi-metal with a spectral weight near momenta {±π/2, ±π/2}, to a large Fermi surface.
	
	Fermi surface reconstruction without symmetry breaking
    Subir Sachdev
    December 30, 2023
======================================================



§ INTRODUCTION

	
	
	Quantum phase transitions involving a change in the volume enclosed by the Fermi surface play a fundamental role in correlated electron compounds. 
	In the cuprates, there is increasing evidence of a phase transition from a low doping pseudogap metal state with small density of fermionic quasiparticles, to a higher doping Fermi liquid (FL) state with a large Fermi surface of electronic quasiparticles <cit.>. In the heavy fermion compounds, much attention has focused on the transitions between metallic states distinguished by whether the Fermi volume counts the localized electronic f moments or not <cit.>. 
	
	Given the strong coupling nature of such transitions, quantum Monte Carlo simulations can offer valuable guides
	to understanding the consequences for experimental observations. As the transitions involve fermions at non-zero density, the sign problem is a strong impediment to simulating large systems. However, progress has been possible in recent years by a judicious choice of microscopic Hamiltonians which are argued to capture the universal properties of the transition but are nevertheless free of the sign problem. Such approaches have focused on density wave ordering transitions <cit.>, where spontaneous translational symmetry breaking accompanies the change in the Fermi volume: consequently both sides of the transition have a Luttinger volume Fermi surface, after the expansion of the unit cell by the density wave ordering has been taken into account. 
	
	Our paper will present Monte Carlo results for quantum phase transitions without symmetry breaking, accompanied by a change in the Fermi surface size from a non-Luttinger volume to a Luttinger volume. 
	The phase with a non-Luttinger volume Fermi surface must necessarily have topological order and emergent gauge degrees of freedom <cit.>.
	There have been a few quantum Monte Carlo studies of fermions coupled to emergent gauge fields <cit.>, and
	our results are based on a generalization of the model of refs. <cit.>, containing a ℤ_2 gauge field coupled to spinful `orthogonal' fermions f_α (spin index α=↑, ↓) carrying a ℤ_2 gauge charge hopping on the square lattice at half-filling. The parameters are chosen so that the f_α experience an average π-flux around each plaquette; consequently when the fluctuations of the ℤ_2 gauge flux are suppressed in the deconfined phase of the ℤ_2 gauge theory, the fermions have the spectrum of massless Dirac fermions at low energies. There are four species of two-component Dirac fermions, arising from the two-fold degeneracies of spin and valley each. In addition to the ℤ_2 gauge charges, the f_α fermions also carry spin and global U(1) charge quantum numbers, and so are identified with the orthogonal fermions of Ref. <cit.>, and the deconfined phase of the ℤ_2 gauge theory is identified as an orthogonal semi-metal (OSM). 
	
	
	It is important to note that the OSM phase preserves all the symmetries of the square lattice Hamiltonian, and it realizes a phase of matter which does not have Luttinger volume Fermi surface. This is compatible with the topological non-perturbative formulation of the Luttinger theorem (LT) <cit.> because ℤ_2 flux has been expelled (about a π-flux background), and the OSM has ℤ_2 topological order <cit.>. 
	There is a Luttinger constraint associated with every unbroken global U(1) symmetry <cit.>, stating that the total volume enclosed by Fermi surfaces of quasiparticles carrying the global charge (along with a phase space factor of 1/(2 π)^d, where d is spatial dimension) must equal the density of the U(1) charge, modulo filled bands. 
	In Oshikawa's argument <cit.>, this constraint is established by placing the system on a torus, and examining the momentum balance upon insertion of one quantum of the global U(1) flux: the Luttinger result follows with the assumption that the only low energy excitations which respond to the flux insertion are the quasiparticles near the Fermi surface. When there is ℤ_2 topological order, a ℤ_2 flux excitation (a `vison') inserted in the cycle of the torus costs negligible energy and can contribute to the momentum balance: consequently, a non-Luttinger volume Fermi surface becomes possible (but is not required) in the presence of topological order <cit.>. In the OSM, the orthogonal fermions f_α carry a global U(1) charge and have a total density of 1: so in the conventional Luttinger approach, there must be 2 Fermi surfaces (one per spin) each enclosing volume (1/2) (2 π)^2. However, with ℤ_2 topological order, the OSM can evade this constraint, and the only zero energy fermionic excitations are at discrete Dirac nodes, and the Fermi surface volume is zero. Earlier numerical studies <cit.> presented indirect evidence for the existence of such a Luttinger-violating OSM phase, and we will present direct evidence here in spectral functions.
	
	Our interest here is in quantum transitions out of the OSM, and in particular, into phases without ℤ_2 topological order. In the previous studies <cit.>, the ℤ_2 confined phases broke either the translational or the global U(1) symmetry. In both cases, there was no requirement for a Fermi surface with a non-zero volume, and all fermionic excitations were gapped once ℤ_2 topological order disappeared. Here, we will extend the previous studies by including an Ising matter (Higgs) field, τ^z, which also carries a ℤ_2 gauge charge. This allows us to define a gauge-invariant local operator with the quantum number of the electron <cit.>:
	
    c_α = τ^z f_α

	The τ^z matter fields do not carry global spin or U(1) charges. With this dynamic Ising matter field present, it is possible to have a ℤ_2-confined phase which does not break any symmetries, and phase transitions which are not associated with broken symmetries. In particular, the Luttinger constraint implies that any ℤ_2 confined phase without broken symmetries must have large Fermi surfaces with volume (1/2) (2 π)^2 for each spin. Our main new results are measurements of the c_α spectral function with evidence for such phases and phase transitions.
	
	One of the unexpected results of our Monte Carlo study is the appearance of an additional “Deconfined FL” phase: see the phase diagram in <ref>. As we turn up the attractive force between the f_α fermions and the Ising matter (Higgs) field τ^z, we find a transition from the OSM to a phase with a large Luttinger-volume Fermi surface of the c_α, but with the ℤ_2 gauge sector remaining deconfined.  Such a phase does not contradict the topological arguments, which do not require (but do allow) a Luttinger-violating Fermi surface in the deconfined state. Only when we also turn up the ℤ_2 gauge fluctuations do we then get a transition to a “Confined FL” phase, which is a conventional Fermi liquid with a large Fermi surface. Within the resolution of our current simulations, we were not able to identify a direct transition from the OSM to the Confined FL, and we have indicated a multicritical point separating them in <ref>.
	
	It is interesting to note that the Deconfined FL phase has some features in common with `fractionalized Fermi liquid' (FL*) phases used in recent work <cit.> to model the pseudogap phase of the cuprates.
	These phases share the presence of excitations with deconfined ℤ_2 gauge charges co-existing with a Fermi surface of gauge-neutral fermions c_α. There is a difference, however, in that the present Deconfined FL state has a large Fermi surface, while the FL* states studied earlier had a small Fermi surface. Both possibilities are allowed by the topological LT. <cit.>.
	
	In passing, we note the recent study of Chen et al. <cit.>, which also examined a ℤ_2 gauge theory coupled to orthogonal fermions, f_α, and an Ising matter field, τ^z. However, their ℤ_2 deconfined phase is different from ours and earlier studies <cit.>: their phase has a large, Luttinger-volume Fermi surface of the f_α fermions which move in a background of zero flux (such zero flux states were also present in the studies of ref. <cit.>). Consequently, their ℤ_2 confinement transition to the Fermi liquid does not involve a change in the Fermi surface volume, and this weakens the connection to finite doping transitions in the cuprate superconductors. 
	
	The rest of the paper is structured as follows. In <ref> we introduce a lattice realization of the Ising-Higgs gauge theory coupled to orthogonal fermions and discuss its global and local symmetries, in <ref> we determine the global phase diagram of our model using  a sign-problem-free QMC simulation. In particular, we study the structure of the Fermi surface and state of the gauge sector in the different phases, and comment on the nature of the numerically observed quantum phase transitions, in <ref> we present a mean field calculation of the physical fermion spectral function in the OSM phase, and lastly in <ref> we summarize our results, discuss relations to experiments and highlight future directions. 
	
	

§ ISING-HIGGS GAUGE THEORY COUPLED TO ORTHOGONAL FERMIONS

	
	
	

 §.§ Lattice Model

	
	
	
	As a concrete microscopic model for orthogonal fermions, we consider the square lattice model depicted in <ref>. The dynamical degrees of freedom are Ising gauge fields, σ^z_b=±1, residing on the square lattice bonds b={r,η}, with r being the lattice site and η=x̂/ŷ, and two types of matter fields: an Ising field τ^z_r=±1 and a spinful orthogonal fermion f_α,r, with α=↑,↓ labeling the spin index. Both matter fields are defined on the lattice sites.  The dynamics is governed by the Hamiltonian, ℋ=ℋ_+ℋ_τ+ℋ_f+ℋ_c comprising the lowest order terms that are invariant under local Ising gauge transformations, as we detail below. 
	
	The first two terms in ℋ correspond to the standard Ising-Higgs gauge theory <cit.>,
	
    ℋ_   =-K∑_□∏_b∈□σ^z_b-g ∑_b σ^x_b  
    ℋ_τ   =-J∑_r,ησ^z_r,ητ^z_r τ^z_r+η -h∑_r τ^x_r.

	In the above equations, the operators σ={σ^x,σ^y,σ^z} and  τ={τ^x,τ^y,τ^z} are the conventional Pauli matrices, acting on the Hilbert spaces of the Ising gauge fields and Ising matter fields, respectively. The Ising magnetic flux term ,Φ_□=∏_b∈□σ_b^z, in ℋ_ equals to the product of the Ising gauge field belonging to the elementary square lattice plaquettes, □. ℋ_τ is a transverse field Ising Hamiltonian for the Ising matter field, where in order to comply with Ising gauge invariance the standard Ising interaction is modified to include an Ising gauge field, σ^z_b, along the corresponding bonds b.
	
	The fermion dynamics is captured by the last two terms in ℋ,
	
    ℋ_f   =-w∑_r,η,ασ^z_r,ηf^†_r,α f_r+η,α+h.c+U∑_r (n^f_r,↑-1/2)(n^f_r,↓-1/2)
    ℋ_c   =-t∑_r,η,ατ^z_rf^†_r,ατ^z_r+ηf_r+η,α+h.c.

	Here, ℋ_f includes a gauge invariant nearest-neighbor hopping of orthogonal fermions and on-site Hubbard interaction between fermion densities, n^f_r,α=f^†_r,αf_r,α. The last term, ℋ_c, defines nearest-neighbor hopping of physical (gauge-neutral) c_r,α fermions as can be readily verified by substituting <ref>. The model is tuned to half-filling (the chemical potential vanishes).
	
	In relation to past works, the model considered here affords a non-trivial generalization of the ones studied previously in refs. <cit.>. In particular, in contrast to prior studies, where the Ising matter fields τ^z were infinitely massive (h→∞ in <ref>), here, varying the transverse field, h, in <ref> controls the excitation gap for τ^z particles. Consequently, the Ising matter fields and subsequently the physical (gauge-neutral) fermion c_α=f_ατ^z are now dynamical degrees of freedom. This important extension provides access to a more generic phase diagram and observables that probe FL physics.
	
	

 §.§ Global and local symmetries

	
	We now turn to discuss the global and local symmetries of our model. ℋ is invariant under a global SU_s(2) symmetry of spin rotations. Furthermore, because we tune to half filling, particle-hole symmetry enlarges the U(1) symmetry, corresponding to fermion number conservation, to form a SU_c(2) pseudo-spin symmetry <cit.>. Physically, the SU_c(2) symmetry generates rotations between the charge density wave and s-wave superconductivity order parameters.
	
	The gauge structure of our model is manifest through the invariance of ℋ under an infinite set of local  gauge transformations generated by the operators G_r=(-1)^n^f_rτ^x_r∏_b∈ +_rσ^x_b. Here, n^f_r=∑_α n^f_r,α and +_r denotes the set of bonds emanating from the site r. Because [ℋ,G_r]=0 for all sites r, the eigenvalue Q_r=±1 of G_r are conserved quantities. Physically, Q_r is identified with the static on-site  background charge assignment. 
	
	To properly define a gauge theory, one must fix the background charge configuration, Q_r. This procedure enforces an Ising variant of Gauss's law G_r=Q_r. Requiring a translationally invariant configuration, two distinct gauge theories may be defined: an even lattice gauge theory with a trivial background (Q_r=1) and an odd lattice gauge theory with a single Ising charge at each site (Q_r=-1). For concreteness, in what follows, we will only consider the case of an odd lattice gauge theory. As explained in ref. <cit.>, at half-filling, the corresponding results for the even sector may be obtained by applying a partial particle-hole transformation acting on one of the spin species <cit.>.
	
	Our model is also invariant under discrete square lattice translations. The operators T̂_x and T̂_y generate translation by a lattice constant along the x and y directions, respectively. When acting on fractionalized excitations, such as the matter fields τ^z_r and f_r,α in our case, translations may be followed by a  gauge transformation. The symmetry operation then forms a projective representation. In the general case, this allows for a richer group structure than the standard (gauge-neutral) linear representation <cit.>.
	
	For the specific case of a  gauge symmetry on a square lattice, a projective implementation of translations is potentially non-trivial. In particular, lattice translation along the x and y directions may either commute or anti-commute <cit.>, namely T_x T_y=± T_y T_x. Physically, the former corresponds to trivial translations, whereas the latter defines a π-flux pattern threading each elementary plaquettes of the square lattice. While for every given choice of a gauge fixing condition, the π-flux lattice inevitably breaks lattice translations, as it leads to doubling of the unit cell, for fractionalized excitations, translational symmetry is restored by applying an Ising gauge transformation <cit.>. This key observation allows for the OSM phase to violate LT without breaking of translational symmetry.
	
	

§ QUANTUM MONTE CARLO SIMULATIONS

	
	

 §.§ Methods

	
	Our model is free of the numerical sign-problem. We can, therefore, elucidate its phase diagram using an unbiased and numerically exact (up to statistical errors) QMC calculations. To control the Trotter discretization errors, we set the imaginary time step to satisfy Δτ≤1/(12|t|), a value for which we found that discretization errors are sufficiently small to obtain convergent results. We explicitly enforce the Ising Gauss law using the methods introduced in ref. <cit.>. Additional details discussing the implementation of the auxiliary-field QMC algorithm and its associated imaginary time path-integral formulation are given in <ref>. Similar results were obtained without imposing the constraint and using the and using the algorithms for lattice fermions (ALF) library  <cit.>.
	
	

 §.§ Observables

	
	To track the evolution of the c electron Fermi surface, we will study the imaginary-time two point Green's function, 𝒢^β_α(k,τ)=-𝒯 [c_k,α(τ)c_k,α^†(0)], where 𝒯 denotes time ordering and the operator c^†_k,α=∑_r f^†_r,ατ^z_re^ik· r creates a c fermion carrying momentum k and spin polarization α. Expectation values are taken with respect to the thermal density matrix, O=1/𝒵[e^-βℋO], with 𝒵=[e^-βℋ] being the thermal partition function at inverse temperature β=1/T. We emphasize that in contrast to the orthogonal f electron, for which, in the absence of a string operator, gauge invariance requires that the two-point Green's function must vanish for all non-equal space-time points, the c electron is a gauge-neutral operator and hence its associated spectral function may be non-trivial.
	
	Determining the Fermi surface structure requires knowledge of real time quantum dynamics. Therefore, some form of analytic continuation of the imaginary time QMC data to real frequency must be carried out. Quite generically, devising a reliable and controlled numerical analytic continuation technique is an outstanding challenge due to the inherent instability of the associated inversion problem <cit.>.
	
	To overcome this difficulty, we employ a commonly used proxy for the low frequency spectral response <cit.>. More explicitly, by computing G^β(k,τ) (the spin index was omitted for brevity) at the largest accessible imaginary time difference τ=β/2, we obtain an estimate for the single particle residue Z. To see how to relate this quantity to real time dynamics, we consider the integral relation 
	
    𝒢^β(k,τ=β/2)=-∫_-∞^∞ d ω A(k,ω)/2cosh(βω),

	where A^β(k,ω)=-1/π𝒢^β(k,iω_m=ω+iδ) is the c fermions spectral function. Since cosh(βω)^-1 tends to unity for βω≪1 and rapidly vanishes in the opposite limit βω≫ 1 , G(k,τ=β/2) amounts to an integral over the spectral function A(k,ω) over a frequency window of order ∼ T. Further assuming a well-behaved spectral response for frequencies ω<T or equivalently no additional low energy excitations, we may use Z̃(k)=β G(k,τ=β/2) as an estimate for the c electron single particle residue Z=A^β(k,ω=0).
	
	To detect the presence of Dirac fermions, we study the finite-size scaling properties of the superfluid stiffness, ρ_s, which captures the long-wavelength transverse current density J_μ(q,ω_m) response to an external electro-magnetic gauge field 𝒜_ν(q,ω_m). Within linear response theory,  J_μ(q,ω_m)=Π_μν(q,ω_m) 𝒜_ν(q,ω_m). Focusing on the ν=μ=x component, the electromagnetic response function equals   <cit.>
		
    Π_xx(q,ω_m)=-[-K_x -J_x(q,ω_m) J_x(-q,-ω_m)].

		Here, K_x is the kinetic energy density associated with x oriented links and the second term is the current-current correlation function.
		In our specific case,
		
    K_x=-∑_r,α w σ^z_r,xf^†_r,α f_r+x̂,α+t τ^z_rf^†_r,ατ^z_r+x̂f_r+x̂,α+h.c.,

		and the current operator at site r along the x direction is given by,
		
    J_r,x=-i[∑_α w σ^z_r,xf^†_r,α f_r+x̂,α+t τ^z_rf^†_r,ατ^z_r+x̂f_r+x̂,α-h.c.].

		With the above definition, we can compute ρ_s as the limit,
		
    ρ_s = lim_q_y → 0ρ_s(q_y) =lim_q_y → 0Π_xx(q_x=0,q_y,ω_m=0)

	
	It is convenient to express the superfluid stiffness in terms of the orbital magnetic susceptibility χ(q) = -ρ_s(q)/q^2. The singularity associated with the Dirac node leads to a diverging magnetic response at low momenta χ(q)=-g_v g_s v_f/(16|q|) <cit.>, with g_v (g_s) being the valley (spin) degeneracy and v_f is the Dirac fermion velocity. Consequently, at low momenta, ρ_s(q_y)∼ q_y ∼ 1/L, so that ρ_s exhibits a slow decay that is inversely proportional to the system size. By contrast, a Fermi liquid admits a finite diamagnetic response at zero momentum (Landau diamagnetism), such that the expected behavior is ρ_s(q_y)∼ q_y^2 ∼ 1/L^2, which vanishes more rapidly than before.
	
	
	To probe a potential instability towards an antiferromagnetic (AFM) order, we monitor the finite-momentum equal-time spin fluctuations χ_S(k)=( ∑_re^ik· r S^z_r)^2, where S^z_r=n^f_↑-n^f_↓ is the f electron spin polarization along the z axis. From χ_S(k), we can compute the staggered magnetization M_AFM=√(χ_S(G_AFM)/L^2), with G_AFM={π,π} being the Bragg vector associated with AFM order. The zero temperature AFM order parameter is obtained by taking the thermodynamic limit Δ_AFM=lim_L→∞lim_β→∞ M_AFM(L,β). In practice, for a given system size, we monitored the convergence of the staggered magnetization toward its zero temperature value. Following that, we extrapolated the finite size data to the infinite system size value, using a polynomial fit in powers of 1/L.
	
	In the presence of dynamical matter fields, determining whether the gauge sector is confined or deconfined is a particularly challenging task due to charge screening. Standard methods, relying on evaluating Wilson loops, no longer sharply distinguish between the two phases. Alternative methods based on extracting the topological contribution to the entanglement entropy <cit.> and the Fredenhagen–Marcu order parameter <cit.> are difficult to reliably scale with system size in fermionic systems. Instead, following refs. <cit.>, we probe the thermodynamic singularity associated with the confinement transition by tracking the Ising flux susceptibility χ_B=∂Φ/∂ K, a quantity that is expected to diverge at the confinement-deconfinement transition, akin to the specific heat singularity in classical phase transitions.
	
	
	
	
	

 §.§ Numerical Results

	
	
	
	
	
	
	
	To render the numerical computation tractable, we must restrict the relatively large parameter space spanned by the set of coupling constants appearing in ℋ. Starting from either the deconfined (g≪ K) or confined phase (g≫ K) our goal is to probe the emergence of physical c fermions and the resulting formation of FL phases in the limit of large hopping amplitude, t. To that end, we numerically map out the phase diagram as a function of the transverse field, g, and physical fermion hopping amplitude, t. Throughout, we will consider a negative Ising flux coupling constant, K<0, for which a π-flux lattice is energetically favorable in deconfined phases.
	
	More concretely, we fix the microscopic parameters w=-1,J=0.1,h=1.0,U=0.1, and K=-1. All energy scales are measured in units of |K|. We note that we have chosen h to be sufficiently large compared to J in order to avoid condensation of Ising matter fields. The resulting two-parameter phase diagram is depicted in <ref>.
	
	We begin our analysis by examining the limiting case t≪ g. In this regime, together with the above choice of microscopic parameters, the Ising matter field τ^z is gapped, and consequently also the physical fermion, c, is expected to be gapped. The low energy physics then involves only orthogonal fermions coupled to a fluctuating Ising gauge field. This physical setting was already studied extensively in previous works <cit.>. 
	
	In the context of our problem, we expect to find a similar structure of quantum phases in the above parameter regime: (i) A confining phase (g
	≫ K), where the orthogonal fermions together with the on-site background static charge (we consider an odd lattice gauge theory) form a localized gauge-neutral bound state, leaving the electronic spin as the only dynamical degree of freedom. Subsequently, quantum fluctuations will generate an effective antiferromagnetic Heisenberg coupling, leading to AFM order at zero temperature. (ii) In the deconfined phase (g
	≪ K), on the other hand, the orthogonal fermions are free and their dispersion is determined by the background flux configuration. For the case of a π-flux lattice, the band structure consists of two gapless and linearly dispersing bands.
	
	To numerically test the above reasoning, in <ref>, we fix t=0.2 and plot the evolution of the flux susceptibility and staggered magnetization as a function of g. Indeed, in agreement with refs. <cit.>, we can identify the aforementioned phases: a deconfined OSM phase for small g, and with an increase in g, we observe a transition towards a confining phase accompanied with AFM order. We use a finite size scaling analysis to estimate the location of both the confinement and AFM symmetry breaking transitions. The flux susceptibility, see <ref>, develops a peak at g_c=0.75(5) that increases with system size and marks the position of the confinement transition. Concomitantly, in <ref>, we find that the AFM order parameter, Δ_AFM begins to rise at g_c=0.75(5).
	
	Due to the increased complexity of the model considered in this work, we found it challenging to reliably estimate universal data associated with the OSM confinement transition, such as critical exponents. Hence, we were unable to make a direct comparison with previous works. Nevertheless, the key signature of the OSM confinement transition, namely the non-trivial co-incidence of confinement and symmetry breaking <cit.> is fully consistent with the numerical data.
	
	
	
	We now turn to address the main inquiry of this study, namely the emergence of low energy gauge-neutral c fermions and their associated spectral signatures. With that goal in mind, in <ref>, we depict our numerical estimate for the momentum resolved c electron residue, Z̃(k), at g=0.55 and for several increasing values of t, beginning from the OSM phase. Remarkably, we find that, in the OSM phase, Z̃(k), comprises four maxima located at momenta k={±π/2,±π/2}. This result is at odds with the conventional LT, which, at half-filling and in the absence of topological order or translational symmetry breaking, predicts a large Fermi surface encompassing half of the Brillouin zone. At large t values, it is energetically favorable for the orthogonal fermion and τ particle to form a gauge-neutral bound state, identified with the c electron, which effectively decouples from the gauge sector. Indeed, with an increase in t, the spectral function continuously evolves into the standard diamond shaped Fermi surface in compliance with LT.
	
	
	
	It is tempting to identify the observed “Fermi pockets” with the nodal points of the OSM. This explanation, however, is incorrect because the orthogonal f fermion is not a gauge invariant object, and in particular the location of the Dirac nodes in momentum space is a gauge-dependent quantity. In <ref>, we provide a simple explanation for this phenomenon using a mean field calculation of the spectral function in the background of a static π-flux configuration.
	
	
	To further track the appearance of c electrons, in Fig. <ref>, we set g=0.55 and study the evolution of the single particle residue evaluated at the anti-nodal point, Z̃(0,π), along the path connecting the OSM and deconfined FL phases as a function of t.  Indeed, we observe that the spectral weight vanishes for t<t_c≈2.0 and continuously rises for t>t_c, signaling the appearance of a large Fermi surface.  
	
	Next, we examine how the emergence of a finite spectral weight for the physical fermions c influences the Dirac orthogonal fermions. To that end, in <ref> we examine the finite-size scaling behavior of the superfluid stiffness by plotting Lρ_s as a function of t for g=0.55. We find that for t<t_c, curves corresponding to different system sizes collapse to a single curve, namely the superfluid stiffness follows a critical scaling ρ_s∼ 1/L. As discussed before, this scaling behavior is characteristic of a low energy gapless Dirac spectrum. Unexpectedly, for t>t_c, a regime in which we previously found a finite quasiparticle weight at the anti-nodal point, we do not observe the expected Fermi-Liquid scaling ρ_s∼ 1/L^2, but rather a behavior consistent with ρ_s∼ 1/L. Although, we can not exclude the possibility that this behavior is due to a finite size crossover, our numerical results indicate the presence of an intermediate phase where orthogonal Dirac fermions and a FL of physical fermions coexist. In Section <ref>, we provide a candidate theoretical description of this scenario.  Lastly, for t ≈ 5.0, the product  L ρ_s vanishes, signaling the absence of low energy orthogonal fermions. The parameter regime t>5 is beyond our current numerical capabilities. 
	
	To further examine the transition, we test for the development of AFM order or confinement in the gauge sector. In <ref>, we  track the evolution of the staggered magnetization as a function of the hopping amplitude, t (with g=0.55, as before). We do not find numerical evidence for a finite AFM order, even when a large Fermi surface is fully developed at large t. Moving to the gauge sector, in <ref>, we probe χ_B  along the same trajectory as above. The flux susceptibility appears to cross the transition smoothly. We can, therefore, conclude that the Fermi surface reconstruction involves neither translational symmetry breaking nor the loss of topological order. Thus, the phase at large t is a deconfined FL.
	
	We remark that due to the perfect nesting condition of the half-filled square lattice, the ground state is expected to exhibits AFM order for arbitrarily small Hubbard interaction. This phenomenon was not observed in our simulations, since in the weak coupling regime U≪ t the magnetization is exponentially small in the coupling constant, and its detection is a notoriously difficult numerical task. 
	
	Next, we examine the path connecting the confined AFM state and the FL phase, by setting g=1.2 and probing the evolution of the spectral function as a function of t. The results of this analysis are shown in <ref>. We observe a featureless flat spectrum deep in the confined AFM phase, as expected due to the absence of fermionic quasiparticles. By contrast, with an increase in t, a large Fermi surface appears. 
	
	To better appreciate the above result, we note that in a confining phase, the low energy spectrum must contain solely gauge-neutral excitations. Indeed, in the confined AFM phase, we can identify these excitations with spin-waves. However, fractionalized orthogonal fermions, carrying U(1) electromagnetic charge, are localized. On the other hand, with an increase in t the attractive force between τ^z and f allows forming a gauge-neutral bound states of c fermions and a metallic state that supports both spin and charge excitations at low energies.
	
	
	We now turn to study the confinement transition along a path connecting the deconfined and confined FL phases. To detect the confinement transition, in <ref>, we plot the flux susceptibility at t=3.0 as a function of g. Indeed, we observe a divergence in χ_B at g_c=0.85(5) marking the location of the confinement transition. The above result demonstrates that our model sustains FL phases in the background of either a confined or a deconfined gauge sector.
	
	Lastly, we summarize our spectral analysis in <ref>, where we plot the spectral weight, Z̃(0,π), evaluated on the “anti-nodal" point as a function of the hopping amplitude t for several values of g. As a starting point, at low t, We consider both the OSM and  AFM phases.  We find that for all g values the spectral weight is small at low t and rises continuously starting from  the energy scale t_c(g). Operationally, we numerically estimate t_c(g) by locating the hopping amplitude for which Z̃(0,π)>0.01. This analysis was used to mark the phase boundaries appearing in <ref>. We remark again, that at half-filling the Fermi-liquid phase is unstable towards the formation of AFM order. Consequently, at strictly zero temperature Z̃(0,π) must vanish for all t and the transition between the AFM and Fermi liquid phases is a smooth cross-over. 
	
	
	
	
	
	
	
	
	

 §.§ Quantum phase transitions

	
	
	We now briefly remark on the theoretical expectations for the quantum phase transitions described above and appearing in Fig. <ref>.
	
		
  * The theory for the transition from the OSM to the AFM was discussed in some detail in ref. <cit.>, and identified as a deconfined critical point with an emergent SO(5) symmetry.
		
  * Away from half-filling, the AFM to Confined FL transition is a conventional symmetry-breaking transition between two confining phases, and is expected to be described by Landau-Ginzburg-Wilson theory combining with damping from Fermi surface excitations, i.e. Hertz-Millis theory <cit.>.
		
  * The transition from the Deconfined FL to the Confined FL phase is a confinement transition without a change in the size of the Fermi surface. It is therefore expected to described by the condensation of an Ising scalar, which can be viewed as representing either the `vison' of the Ising gauge theory or the Ising matter field τ^z <cit.>. The large Fermi surface of electron-like quasiparticles will damp the quasiparticles, but this damping is much weaker than that in Hertz-Millis theory: the resulting field theory was described in refs. <cit.>. We note that this field theory also applies to the confinement transition in ref. <cit.>, where the damping is due to a large Fermi surface of orthogonal fermions, and there is no change in the size of the Fermi surface across the transition.
		
  * We do not have a theory for a direct transition from the OSM to the Deconfined FL with a large Fermi surface transition. Indeed,  it may well be that this transition occurs via an intermediate phase, with gapless excitations of both electrons and orthogonal fermions <cit.>. The presence of a regime within the “Deconfined FL” region of the phase diagram where we appear to observe co-existence of a c Fermi surface (as measured by the c spectral density) with a ρ_s which vanishes as ∼ 1/L is evidence in support of such an intermediate phase. The Luttinger constraint requires that the orthogonal fermions which are absorbed into the c Fermi surface must be accounted for by those forming the Dirac nodes. One way to preserve the Dirac nodes in an intermediate phase, and also maintain consistency with particle-hole symmetry, is to form electron-like and hole-like pockets of an equal area of the c fermions; the c hole pockets would then be in ancillary trivial insulator, similar to refs. <cit.>. Our resolution is not sharp enough to resolve such an intricate Fermi surface evolution, which is likely present in the intermediate phase, and we leave its study to future work.
	
	
	

§ MEAN FIELD CALCULATION OF THE PHYSICAL FERMION SPECTRAL FUNCTION IN THE OSM PHASE.

	
	
	The numerical observation of a finite spectral weight centered about the four nodal points k={±π/2,±π/2} is surprising and requires further analytic understanding. To that end, in this section, we will present a simple and intuitive explanation using a mean-field calculation of the spectral function of c fermions in the OSM phase. In our calculation, we consider a static background π-flux configuration and neglect gauge field fluctuations. This approximation is justified deep in the OSM phase, where such fluctuations are small due to the finite vison gap. In this setting, we can model the dynamics of the f electron ( τ^z field) by a free fermion (scalar field) hopping in the background of a π-flux lattice.
	
	To make progress, we must choose a concrete a gauge fixing condition for the Ising gauge field. We take an Ising variant of Landau's gauge, σ_r,x̂=(-1)^r_y and σ_r,ŷ=1, see <ref>. We emphasize that while the gauge fixing procedure inevitably breaks both lattice translations and π/2 rotations, expectation values of physical (gauge-neutral) observable must preserve these symmetries. More concretely, we take the effective Hamiltonians,
	
    ℋ^MF_f   =-t_f∑_r,η t_r,η f_r^† f_r+η +h.c.
    ℋ^MF_ϕ   =∑_rπ_r^2/2m +m ω^2/2( ∑_r Δϕ_r^2 +∑_r,η(ϕ_r-t_r,ηϕ_r+η)^2),

	
	Here, ϕ_r is a scalar real field and π_r is its canonically conjugate momentum, [ϕ_r,π_r']=iδ _r,r'. With the above gauge choice the hopping amplitudes are set to t_r,η=((-1)^r_yδ_η,x̂+δ_η,ŷ). The scalar field Hamiltonian is parameterized by the inertial mass m, oscillation frequency ω, and Δ, which allows controlling the single particle excitation gap for ϕ particles.
	
	We tune Δ to work in a regime where, on the one hand, the scalar field has a finite gap to avoid condensation, but, on the other hand, it is sufficiently small to render the physical fermion gap small. The lowest order diagram <cit.> contributing to 𝒢^β_c(k,iω_m) is the bubble diagram shown in  <ref>. Evaluating the diagram boils down to a convolution of the Ising matter field propagator, 𝒢^β_ϕ(k,iω_m), and the orthogonal fermion propagator, 𝒢^β_f(k,iω_m). Both propagators can be readily computed by diagonalizing the quadratic Hamiltonians in  <ref>. Further details of this calculation are given in <ref>.
	
	In <ref> we depict the zero-frequency spectral function A(k,ω=0)=-1/π[𝒢^β_c(k,iω_m=i0^+)], evaluated for the microscopic parameters t=m=ω=1 and Δ=-1.1, and inverse temperature β=1/32.
	Remarkably, our simplified model displays a finite spectral weight located at the nodal points k={±π/2,±π/2}, in agreement with the exact QMC calculation. As a non-trivial check of our computation, we observe that unlike the gauge fixed hopping Hamiltonians of <ref>, the physical spectral function respects the full square lattice C_4 point group symmetry. We note that at strictly zero temperature, the spectral weight must vanish due to the non-zero Ising matter field gap.
	
	An intuitive understanding of this result may be derived by examining the Dirac band structure in the π-flux phase, shown as a contour plot in  <ref>. At low temperatures, the orthogonal fermions will occupy all k-space modes up to the Dirac point at {±π/2,π/2}. However, the Ising matter field, τ^z follows Bose statistics and hence will concentrate at the band minimum, {0/π,0}. As a result, the integral over the internal momentum in <ref> will be appreciable only at momentum transfer k={±π/2,±π/2} (black arrow in <ref> connecting the f and τ^z particle). In other words, the momentum-space splitting of the two flavors of matter fields is responsible for the unconventional spectral response of the OSM phase.
	
	
	

§ DISCUSSION AND SUMMARY

	
	
	We have studied a lattice model of orthogonal-fermions coupled to an Ising-Higgs gauge theory. The absence of the sign-problem enabled us to determine its global phase diagram and explore related phase transitions using a numerically exact QMC simulation. A key ingredient of our study, which non-trivially distinguishes it from previous works, is the introduction of an Ising matter field that together with the orthogonal fermion may form a physical fermion. This crucial feature of our model enabled access to the study of FL phases in the presence of  topological order.
	
	Notably, our model hosts both non-FL quantum states that violate LT due to  topological order, and also LT preserving FL states, where the gauge sector is either confined or deconfined. On tuning of microscopic parameters, we are able to cross quantum critical points that separate these phases, some of which appear to be continuous.
	
	It is interesting to make a connection, even if suggestive, between our numerical results and experimental signatures of Fermi surface reconstruction observed in cuprate materials. In particular, the OSM phase shares several properties with the pseudo-gap phase: a strong depletion in the density of states at the anti-nodal points, and concentration of finite fermionic spectral weight at the nodal points k={±π/2,±/2}. Most importantly, there is experimental evidence that these phenomena occur in the absence of translational symmetry breaking. We also described phase transitions involving the appearance of a large Fermi surface from such a state, and this has connections to phenomena in the cuprates near optimal doping <cit.>.
	
	
	Our study has left some open questions. In particular, it would be interesting to develop a field theory description to the transition between the OSM and deconfined FL, which numerically appears to be continuous. Such a description will have to address the unusual finite-size scaling of the superfluid stiffness in the deconfined FL phase.  From the numerical perspective, simulations on larger lattices are key in resolving the properties of this transition. In addition, it would be interesting to study the fate of the OSM phase and neighboring phases away from half-filling. This can be achieved, using a sign-problem free QMC, at least for the even lattice gauge theory case. We leave these interesting questions for future studies. 
	
	
	A note added: During the final stages of this work, we have become aware of a related work by Chen et al. <cit.> studying a lattice model of orthogonal fermions in a parameter regime complementary to our work, which we described briefly in <ref> and <ref>.
	
	
		We thank Ashvin Vishwanath, Aavishkar A. Patel, and Amit Keren for useful discussions. S.G. acknowledges
		support from the Israel Science Foundation, Grant No. 1686/18.
		S.S. was supported by the National Science Foundation under Grant No. DMR-1664842. S.S. was also supported by the Simons Collaboration on Ultra-Quantum Matter, which is a grant from the Simons Foundation (651440, S.S.)
		F.F.A. thanks the DFG collaborative research centre SFB1170 ToCoTronics (project C01) for financial support as well as the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter – ct.qmat (EXC 2147, project-id 390858490).
		This work was partially performed at the Kavli Institute for Theoretical Physics (NSF grant PHY-1748958). FFA gratefully acknowledges the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SUPERMUC-NG at the Leibniz Supercomputing Centre (www.lrz.de). This research used the Lawrencium computational cluster resource provided by the IT Division at the Lawrence Berkeley National Laboratory (Supported by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231) and the Intel Labs Academic Compute Environment.
	
	
	

§ PATH INTEGRAL FORMULATION

	
	The partition function, 𝒵(β) at inverse temperature β is given by,
	
    𝒵(β)=Tr[P̂e^-βℋ]

	where the projection operator P̂=∏_iP_i with,
	
    P̂_i=1/2(1+(-1)^n^f_iτ^x_i∏_+σ_ij^x)
 
	We rewrite the projection operator using a discrete Lagrange multiplier as,
	
    P̂_i=∑_λ_i±1P̂_i,λ=∑_λ_i=±1e^iπ/2 (1-λ_i)(∑_+(1-σ^x_i,j/2)+1-τ_i^x/2 +n^f_i)

	
	Next, we use Trotter decomposition to write e^-βℋ=∏_m=0^M-1 e^-ϵℋ with ϵ=β/M and insert a resolution of identities in the σ^z and τ^z basis, 1=∑|τ^z,σ^z⟩⟨τ^z,σ^z| leading to,
	
    𝒵(β)=∑_λ,τ^z,σ^zTr_f ⟨τ^z_0,σ^z_0|∏_iP̂_i,λ e^-ϵℋ|τ^z_M-1,σ^z_M-1⟩…⟨τ^z_1,σ^z_1|  e^-ϵℋ|τ^z_0,σ^z_0⟩

	
	We first focus on the last time step, which contains P̂_i,λ. In the Ising sector, the only off-diagonal (in the σ^z,τ^z basis ) terms are the transverse field and the constraint. Focusing on a specific site i, for the  matter field we obtain the matrix element,
	
    ⟨τ^z_i,0|  e^iπ/2 (1-λ_i)(1-τ_i^x/2)+ϵ h τ_i^x|τ^z_i,M-1⟩   =∑_τ^x_i=±1⟨τ^z_i,0|  e^iπ/2 (1-λ_i)(1-τ_i^x/2)+ϵ h τ_i^x|τ^x_i⟩⟨τ^x_i||τ^z_i,M-1⟩
       =1/2∑_τ^x_i=±1 e^ϵ h τ_i^x e^iπ(1-τ^x_i/2)(1-λ_i/2+1-τ^z_i,0/2+1-τ^z_i,M-1/2)
       =1/2(e^ϵ h+e^-ϵ he^iπ(1-λ_i/2+1-τ^z_i,0/2+1-τ^z_i,M-1/2))
       =cosh(ϵ h)    λ_iτ^z_i,0τ^z_i,M-1=1 
    sinh(ϵ h)    λ_iτ^z_i,0τ^z_i,M-1=-1

	The effective Boltzmann weight is then,
	
    W (λ_i,τ^z_i,0 ,τ^z_i,M-1)∝ e^γτ^z_i,0λ_i τ^z_i,M-1

	where γ=-1/2log(tanh(ϵ h)).
	Physically, the above action corresponds to the gauge invariant Ising interaction along the temporal direction.
	Importantly, we must take h>0 in order to avoid a sign problem.
	
	Similarly to ref.<cit.>, the constraint term associated with the Ising gauge field leads to a spatio-temporal plaquette term in the 3D Ising gauge theory, and the f fermions Green's function is modified by the introduction of a diagonal matrix P[λ_i] with diagonal elements P_ii=λ_i.
	
	
	
	
	

§ DETAILS OF THE MEAN FIELD CALCULATION

	
	
	To evaluate the bubble diagram in <ref>, we first express the mean field Hamiltonians (<ref>) in momentum space defined on a reduced Brillouin zone (0<k_x<2π,0<k_y<π), as imposed by the gauge fixing choice in <ref>. Explicitly, for the fermionic part, substituting f_r=1/√(N)∑_ke^ikrf_k (N being the number of lattice sites) gives
	
	
    ℋ^MF_f   =-t_f∑_k,k' f^†_kf_k'(2cos(k_y)δ_k,k'+2cos(k_x)δ_k,k'+πk̂_y)
       = -t_f∑_0<k_x<2π,0<k_y<π[ f^†_k,0 f^†_k,π ][  2cosk_y  2cosk_x;  2cosk_x -2cosk_y ][ f_k,0; f_k,π ]

	where f^†_k,0=f^†_k and f^†_k,π=f^†_k+πk̂_y. Diagonalizing the above Hamiltonian, we obtain the fermionic eigen spectrum ϵ_±(k) and eigen-modes f_α(k)=V_α,γ(k)f_γ(k), where γ=± and α=0/π and V_α,γ is the diagonalizng matrix.
	
	A similar analysis is carried out in order to diagonalize the scalar field ϕ mean field Hamiltonian. Writing <ref>, in momentum space gives, (ϕ_r/π_r=1/√(N)∑_kϕ_k/π_k e^ikr)
	
    ℋ_ϕ   =∑_kπ_k π_-k/2 + mω^2/2( ∑_kΔϕ_k ϕ_-k+∑_k,k'ϕ_kϕ_k'( 4 - 2 cos(k_y) δ_k,-k'-2cos(k_x)δ_k,-k'+πk̂_y)) 
       =∑_0<k_x<2π,0<k_y<π∑_απ_k,απ_-k,α/2 + mω^2/2∑_α,α'ϕ_k,αK_α,α'(k)ϕ_-k,α',

	where the spring constant matrix equals,
	
    K_α,α'(k)=[ Δ +4-2 cos(k_y)      -2cos(k_x);      -2cos(k_x) Δ +4+2 cos(k_y) ].

	
	By diagonalizaing K_α,α' (the mass matrix is already diagonal) we obtain the eigen-frequencies ω_κ(k) and the normal modes ϕ_k,α=U_α,κ(k)ϕ_k,κ. As before, κ=±, and α=0,π.
	With the above definitions, our final expression of the bubble diagram amplitude reads,
	
    𝒢_c(k,iω_m)    = ∑_q,q',ν_m,γ,κδ_p(q),p(q')δ_p(k-q),p(k-q')V_α(q),γ(p(q)) 
       × V_α(q'),γ(p(q')) U_α(k-q),κ(p(k-q)) U_α(k-q'),κ(p(k-q')) 
       ×1/iν_m-ϵ_γ(p(q))1/((ν_m-ω_m)^2+ω^2_κ(p(k-q))

	
	Here, the reduced momentum p(q) is defined as p(q)={q_x,q_y π} and the momentum index α(q) equals 0 (π) if q_y∈[0,π] (q_y∈[π,2π]), as before. We note that the resulting momentum integration is restricted to the two cases q=q' or q=q'+πq̂_y.
	
	
	The Matsubara sum can be evaluated analytically,
	
    ∑_ν_m1/iν_m-ϵ_γ1/(ν_m-ω_m)^2+ω^2_κ=- βω _κtanh(βϵ _γ/2)+ β(ϵ_γ-i ω_m) (βω_κ/2)/2 ω _κ(ω_κ^2- (ϵ _γ-i ω_m)^2)

	
	and we are left with the momentum integration that is computed numerically on a discretized momentum grid.
	
	
	
