Emergence of Gauss' Law in a $Z_2$ Lattice Gauge Theory

We explore a $Z_2$ Hamiltonian lattice gauge theory in one spatial dimension with a coupling $h$, without imposing any Gauss' law constraint. We show that in our model $h=0$ is a free deconfined quantum critical point containing massless staggered fermions where all Gauss' law sectors are equivalent. The coupling $h$ is a relevant perturbation of this critical point and fermions become massive due to confinement and chiral symmetry breaking. To study the emergent Gauss' law sectors at low temperatures in this massive phase we use a quantum Monte Carlo method that samples configurations of the partition function written in a basis in which local conserved charges are diagonal. We find that two Gauss' law sectors, related by particle-hole symmetry, emerge naturally. When the system is doped with an extra particle, many more Gauss's law sectors related by translation invariance emerge. Using results in the range $0.01<h \leq 0.15 $ we find that three different mass scales of the model behave like $h^{p}$ where $p \approx 0.579$.


Introduction
The possibility of using quantum computers to understand quantum field theories has become an exciting field of research lately [1]. One of the challenges for applications to nuclear and particle physics is our ability to formulate all local quantum fields on a finite dimensional Hilbert space [2]. While there are many ways to accomplish this [3,4], the low energy physics that emerges can depend on the details of the Hilbert space formulation [5]. Simple quantum field theories are currently being formulated so that they can be studied using a quantum computer [6,7,8,9,10,11,12,13]. One of the main long-term challenges is to be able to formulate and study strongly coupled gauge theories like QCD [14]. One of the goals of our work is to understand the physics of a simple qubit-based quantum field theory with similarities to QCD, but simple enough to be simulated on a quantum computer in the near future.
Since the Hamiltonian of a gauge theory is invariant under local symmetry transformations, the Hilbert space of states can be divided into sectors labeled by local conserved gauge charges [15]. Under time evolution these sectors do not mix with each other and each sector satisfies its own Gauss' law. For example, if the local conserved gauge charge density is chosen to be ρ G (r) in quantum Electrodynamics, the Gauss' law will take the form ∇ ∇ ∇ · E(r, t) − ρ(r, t) = ρ G (r). Using Maxwell's equations we learn that the sector of the quantum Hilbert space with ρ G (r) = 0 is the physical one. Imposing this local constraint on the allowed space of states is an important step in the formulation of the quantum gauge theory.
Since imposing the Gauss' law constraint can be difficult when formulating a gauge theory on a quantum computer [9], in this work we explore if such constraints can emerge naturally at low temperatures without imposing them. We study this question using a simple one dimensional Z 2 lattice gauge theory. Such gauge theories in higher dimensions are believed to describe frustrated quantum magnets and spin liquid phases of materials [16,17]. Our model contains massless fermions interacting with a Z 2 lattice gauge field, whose Hamiltonian is given by Here c † j and c j create and annihilate fermions on the sites j = 0, 1.., L − 1 of a periodic lattice and σ x j , σ z j are the Pauli matrices that represent the Z 2 gauge fields associated to links connecting the sites j and j + 1. We assume L to be even to preserve particle hole symmetry.
It is easy to verify that all local charge operators Q j = σ z j−1 σ z j (−1) nj , where n j = c † j c j , commute with each other and with the Hamiltonian. The set of their simultaneous eigenvalues {Q j = ±1} label one of 2 L possible Gauss' law sectors. We label each sector with a unique number from 0 to 2 L − 1 using the relation and compute the probability distribution P (Q) as a function of temperature. The emergent sectors Q e are those which have a non-zero probability P (Q e ) at zero temperature. In the next section we will argue that when h = 0 our lattice model describes a deconfined quantum critical point where all Gauss' law sectors are equivalent. In later sections, using Monte Carlo calculations, we will show that when h = 0 two Gauss' law sectors related by particle-hole symmetry emerge, fermions are confined and acquire a mass due to chiral symmetry breaking. These properties make the model similar to QCD. The emergence of Gauss' law sectors when h = 0 was discussed recently and studied using auxiliary field Monte Carlo (AFMC) methods in two spatial dimensions [18,19]. Unfortunately, it is not possible to compute P (Q) using AFQMC calculations due to sign problems. Since our studies are restricted to one spatial dimension, we can work in the basis in which fermions are represented by their occupation numbers and the Z 2 electric field operators (σ z j ) are diagonal, without encountering sign problems. Since every configuration is naturally in a fixed Gauss' law sector with a well defined set of charges {Q j }, we can easily compute the probability distribution P (Q) using our Monte Carlo method. We can also address the effects of doping the system away from the particle-hole symmetric situation. We show that L different Gauss' law sectors related by translation symmetry emerge if we dope the system with one extra particle.
Our calculations are performed in a discrete time formulation of the path integral, where we divide the inverse temperature β into equal slices of width ε = 0.1 [20]. An illustration of the worldline configuration is shown in Fig. 1. While the algorithm is straight forward [21,22,23], more details about how we overcome some issues that arise in a gauge theory can be found in the supplementary material.

Deconfined Quantum Critical Point
In order to understand our model Eq. (1) better, let us first set h = 0 and ignore the gauge fields in the fermion hopping term. If we then add a staggered fermion mass term H (m) j = (−1) j+nj to the model it is easy to verify that the resulting free fermion Hamiltonian, describes relativistic staggered fermions with mass m [24]. In this sense our model describes the physics of massless staggered fermions coupled to Z 2 lattice gauge fields. In fact we can solve our model exactly when h = 0. To see this, let us define a new set of fermion annihilation operators through the relations The new fermion creation operators will naturally be the Hermitian conjugates and n j = f † j f j = c † j c j . It is easy to verify that these new set of fermion operators non only satisfy the usual anti-commutation relations but also commute with the local gauge charges Q j = σ z j−1 σ z j (−1) nj defined previously. However, it is important to remember that Q j 's satisfy the constraint, where N f is the total fermion number. This shows that a choice of the Gauss' law sector {Q j ± 1}, also constrains the fermionic Hilbert space. Let us also define two new operators in the gauge field sector: the Wilson loop operator W L−1 = σ x 0 σ x 1 ...σ x L−1 , and its conjugate E L−1 = σ z L−1 . Although f j , f † j depend on gauge fields they still commute with W L−1 and E L−1 .
ε Figure 1: Illustration of a worldline configuration of fermions (circles) and gauge fields (dark (σ z j = +1) and light (σ z j = −1) ovals). Transfer matrix elements e εH are divided into fermion hopping (dark plaquettes) and the coupling h (light plaquettes). Our time discretization error is chosen to be ε = 0.1. Note that we can write σ z n k is a non-local fermion operator. Hence, we can rewrite Eq. (1) as where for convenience we define W j = 1 for j = L − 1. Since all fermion operators commute with all W j , Q j and E L−1 , this form of the Hamiltonian makes it easy to analyze the physics in each Gauss' law sector. We can work in a basis where all {Q j } are diagonal and the constraint Eq. (4) is satisfied. When h = 0 we observe that the fermions are almost free and massless, except for the coupling to W L−1 . In a basis where the Wilson loop is diagonal, the choice of W L−1 = ±1 only effects the boundary condition of the free fermion problem and disappears in the thermodynamic limit. Thus, we find that all Gauss' law sectors are equivalent and describe free massless staggered fermions. This implies that h = 0 is a deconfined quantum critical point in every Gauss' law sector.

Emergence of Gauss' Law
When h = 0 different Gauss' law sectors have different energies and as far as we know our model cannot be solved exactly. In our work we use Monte Carlo methods to explore the physics. We first study the emergence of Gauss' law at low temperatures when L = 12 at h = 0.05 by computing the probability distribution P (Q) for Q = 0, 1, ..., 4095 at various inverse temperatures β. In Fig. 2 we plot P (Q) as a function Q at β = 5, 25, 50, 100. As the temperature reduces, we observe two Gauss' law sectors emerge with Q e = 1365 (i.e., {Q j = (−1) j }) and its particle-hole symmetric partner Q e = 2730 (i.e., {Q j = (−1) j+1 }). In these emergent sectors the fermion number is found to satisfy the constraint Eq. (4) and is consistent with particle-hole symmetry i.e., N f = L/2 = 6. These half-filled pattern, {Q j = (−1) j } and {Q j = (−1) j+1 } with N f = L/2 continue to be the emergent sectors even at larger lattice sizes, as predicted in the previous work [18,19].
We then understand how P (Q e ) depends on the temperature and what is the role of the lattice size. In Fig. 3 we plot the P (Q e ) for each of the two emergent sectors as a function of L for different values of β at h = 0.05. Remarkably, although the number of Gauss' law sectors increase exponentially as 2 L , P (Q e ) saturates to the thermodynamic value for L > 24. This implies that in the thermodynamic limit only a few sectors are important at low temperatures. On the other hand very low temperatures are required to project out all the subdominant sectors completely. For example   at h = 0.05 we need β = 250 for P (Q e ) ≈ 0.5. It is not surprising that this temperature is dependent on h since at h = 0 all sectors are equally important. In Fig. 4 we plot P (Q e ) as a function of L at various combinations of (h, β). Note that at h = 0.01, even β = 750 is still not sufficient to project out all the sub-dominant Gauss' law sectors.
Changing N f must change Q e as expected from Eq. (4). To study this, we add a chemical potential term µ to the Hamiltonian and increase it to dope the system with one additional fermion. At h = 0.05, L = 12 and β = 100 we observe that when µ increases from 0 to 0.2 the fermion number increases from N f = 6 to N f = 7 (for a plot of the µ dependence we refer the reader to the supplementary material). When N f = 7 (at µ = 0.2) we observe that twelve Gauss' law sectors given by Q e = 3434, 2773, 1451, 2902, 1709, 3418, 2741, 1387, 2774, 1453, 2906, 1717 emerge at low temperatures. In Fig. 5 we show P (Q) close to these emergent sectors. These sectors are consistent with translation symmetry and satisfy the constraint Eq. (4). The presence of one extra fermion creates two defects in the half filled Gauss' law pattern, that are maximally separated (a pictorial representation of the emergent sectors can be found in the supplementary material).

Chiral Symmetry Breaking
In Section 2 we argued that the local operator H (m) j is a fermion mass term in our model. Thus, the ground state can be viewed as the chiral order parameter for the chiral symmetry that preserves the masslessness of fermions. Note that H m j is odd under both particle-hole transformations and translation by one lattice spacing, both of which are important to maintain the masslessness of fermions. Hence we will refer to both these symmetries together as the chiral symmetry in our model. While our model is invariant under chiral symmetry for all values of h, this symmetry is explicitly broken in a fixed Gauss law sector. For example in each of the two Gauss law sectors {Q j = (−1) j } and {Q j = (−1) j+1 } that emerge at low temperatures when h = 0, the chiral symmetry is explicitly broken and fermions become massive. This implies that fermions will become massive when h = 0. When h = 0 an extra symmetry related to the transformation σ z → −σ z on alternate bonds can be used to redefine chiral symmetry in each Gauss' law sector. This helps to keep fermions massless at the deconfined quantum critical point, but not away from it.
Since we do not impose any Gauss' law constraint, chiral symmetry is never explicitly broken and the fermion mass generation at h = 0 must be viewed as arising due to spontaneous chiral symmetry breaking that can be observed   Table 1. through a non-zero value of φ. In our Monte Carlo calculation, it can be extracted using the chiral susceptibility which is expected to scale as χ = φ 2 L + A in the symmetry broken phase. In Fig. 6 we show our data for χ at various values of (h, β). The solid lines are fits to the expected form in the broken phase and the extracted values of φ and A are tabulated in Table 1. Thus we see that when h = 0, our model describes a chirally broken massive fermion phase. Since the energy of the gauge string connecting two fermions excited over the ground state in a given Gauss' law sector will grow with h, fermions remain confined. Thus at non-zero couplings our model describes massive confined fermions like in QCD.

Critical Exponent and Scaling
The deconfined quantum critical point at h = 0 can be used to define a continuum limit of our massive Z 2 lattice gauge theory [25]. The effective action of the two dimensional continuum quantum field theory (QFT) that emerges will be described a theory of massless fermions with a relevant coupling at leading order of the form  As h vanishes we expect the chiral order parameter to vanish as φ ∼ h p , where p is a critical exponent (anomalous dimension) that depends on the dimensions of φ and O h at the critical point.
The simplest possibility is that O h is just a fermion bilinear mass term. This is consistent with the fact that chiral symmetry is broken in each fixed Gauss' law sector that emerges when h = 0. Since the chiral order parameter depends linearly on the fermion mass, this implies that p = 1. Of course O h will also contain four-fermion couplings, but these would be less relevant and can be ignored in the current discussion. On the other hand, since the lattice Hamiltonian Eq. (5) contains a long range interaction F j , we cannot rule out the possibility that O h contains a term that is more relevant than a mass term. If this is true then p < 1 is also possible.
We can study this question by fitting the values of φ in Table 1 to the form φ ∼ h p . We find an excellent fit for p = 0.578(2) if we drop the h = 0.01 point, which is justified since P (Q e ) has not yet saturated to 0.5 at (h = 0.01, β = 750) (see Fig. 4). Although it is possible that the range of h we have used in the fit is influencing the value of p, the fact that it is considerably smaller than 1 suggests that the long range part of F j is indeed playing an interesting role.
At the free massless fermion fixed point in two dimensions, the chiral order parameter φ has the canonical dimensions of a mass since it is a fermion bilinear. Assuming this argument extends to our deconfined quantum critical point, φ ∼ h p implies that other quantities with the dimension of mass induced by h will also scale similarly. We have verified this scaling prediction using two other quantities with mass dimensions. We first consider the term h σ z , which is just the average of the Hamiltonian density in one dimension and has dimensions of a square of a mass scale. Thus we can define M g = h σ z which has dimensions of a mass, and according to our above prediction we expect M g ∼ h p . In Table 1 we show our Monte Carlo results for σ z at various values of (h, β). If we use this data to extract M g and fit it to the form M g ∼ h p we get p = 0.577 (2), in excellent agreement with our result above.
We can also compute a winding mass M w , obtained from the exponential decay of the spatial winding number susceptibilty W 2 with spatial size L. Here W is the spatial winding of fermion worldlines of our Monte Carlo configurations. We expect W 2 /β = B exp(−M w L) at low temperatures. In Fig. 7 we show our data for W 2 /β as a function of L for four different values of (h, β). The solid lines are fits of the data to the expected form. The values of M w and B extracted from these fits are tabulated in Table 1. Due to an exponentially small signal, we are unable to extract M w reliably for h > 0.1. Fitting the data in the table to M w ∼ h p we get p = 0.589 (8), again in excellent agreement with the previous two results. The data for φ, M g and M w as functions of h are shown in Fig. 8. We find that Ch p , with p ≈ 0.579 describes all the three mass scales in the region shown.

Conclusions
In this work we have studied a simple Z 2 lattice gauge theory where all local degrees of freedom can be represented by single qubits, so that it can be easily explored on a quantum computer. Our model has an interesting deconfined quantum critical point, which when perturbed by a relavant coupling h leads to a massive quantum field theory similar to QCD. The coupling h seems to be more relevant than a simple fermion mass term. Although we did not impose any Gauss' law in our model, it emerges naturally at low temperatures in the massive phase. Doping the system changes the emergent Gauss' law sectors. Extensions of our work to higher dimensions with bosonic matter fields is easy and would also be interesting especially given the richness of phase diagrams of such theories [26].

Supplementary Material Appendix A. Worldline Monte Carlo Algorithm
Here we briefly explain our Monte Carlo algorithm which is a simple extension of well known algorithms developed earlier to update worldline configurations [21,22,23]. In our algorithm we begin with the partition function is essentially our model with an additional chemical potential term. By tuning the value of µ we can change the particle number N f in the ground state. This helps us understand how emergent Gauss' law sectors change due to doping.
In order to express the partition function as a sum over weights of worldline configurations, we first divide β into L T equal slices of width ε and then express the small imaginary time evolution operator as e −εHµ = e −εHo e −εHe , where and H o is a similar operator with the sum over j performed over odd sites. Note that both these operators contain a sum over terms that act on bonds between neighboring sites that commute with each other. Hence one can easily compute the matrix elements of e −εHo and e −εHe as a product of weights on local two dimensional plaquettes as shown in the left figure of Fig. A.10. We distinguish two types of plaquettes: fermion plaquettes (shown with darker shade) and gauge plaquettes (shown with lighter shade). The weight of a fermion plaquette is obtained from the transfer matrix elements where |n j , n j+1 , s j are eigenstates of c † j c j , c † j+1 c j+1 and σ z j . Note that these matrix elements have non-zero entries in both diagonal and off-diagonal terms. On the other hand the weight of the gauge plaquette is based on the matrix elements n j n j+1 s j | exp εhσ z j |n j , n j+1 s j .
that are always diagonal. For this reason the values s j across a gauge plaquette are constrained to be the same. As shown in Fig. A.10 the fermion and gauge plaquettes occur alternately on each time slice. During a worm update, we pick a site at random and propose to place a creation operator on that site (worm tail) and place an annihilation operator on one of the other three sites associated with a fermion plaquette (worm head). If this move is accepted according to detailed balance, we accept it as a new configuration. If the configuration is accepted we pick a neighboring fermion plaquette connected to the site containing the annihilation operator (that can be either forward or backward in time) and perform the same operation of adding a creation and annihilation operator. However, now the creation operator is added to the site connected to the previous plaquette and the annihilation operator is chosen to be one of the remaining three sites. If the proposal is accepted the worm head is moved to the new site and the update proceeds. If the proposal is rejected we go back to the previous plaquette and again the update proceeds. Note that these configurations with a worm head and a tail are not part of the partition function since they contain an extra creation and annihilation operator separated by some distance. We refer to these configurations as being in the worm sector. The worm update may only end when the worm head meets the worm tail and we get an acceptable configuration in the partition function sector. The above steps in constructing the worm update are standard, except that in a gauge theory new complications arise and need to be addressed. The local gauge constraints will be violated at each step unless the gauge links are also simultaneously updated. However, it is impossible to remain in the same Gauss' law sector if only a single fermion is created on a time slice. This is because every fermion creation changes the Gauss' law sector. However, since we allow all Gauss' law sectors to be sampled this is not a problem for our study. On the other hand we do still have to update the gauge links associated to a fermion plaquette whenever the worm head moves to the neighboring site. This is because such a move introduces or elimates a fermion hop and hence will add or remove one σ x j operator between sites j and j + 1. But addition of this operator at one time slice has a non-local effect on all time slices in time. Further the trace constraint of the partition function will be violated.
For this reason in our update we do not impose the trace constraint in the worm sector. We allow the gauge links across a fixed temporal boundary to fluctuate freely. This temporal boundary is chosen to be on the time slice of the tail site for the worm update. A partial worm update along with the temporal boundary chosen is illustrated on the right in Fig. A.10. This then allows us to perform a local worm update in the fermion sector as discussed above, but such an update always includes the possibility of an extra update of all the gauge links associated with σ z j -which is non-local in time-when the fermion plaquette between sites j and j + 1 is updated. This makes our worm update non-local in time but local in space.

Appendix B. Observables and Tests
We have measured four observables in our Monte Carlo calculations. These are

Winding Number Susceptibility
where H θ µ has an extra factor e iθ or e −iθ when a fermion hops across the boundary depending on the direction of the hop. 4. Average fermion number We have tested the algorithm by computing these observables exactly on a L = 4 lattice. In Table B.2 we compare our exact answers with those obtained using the algorithm for ε = 0.01 (referred to as MC1) and for ε = 0.1 (referred to as MC2). Note that there is not much difference between the two results and both of them match well with the exact calculations. For this reason in our work we have fixed ε = 0.1 in our calculations.

Appendix C. Summary of Additional Plots
In order to document some of our results in a more complete form, we provide several clarifying figures.  The values of σ z we use in the main paper are estimated from this data.