{'fraction_non_alphanumeric': 0.036103812638224674, 'fraction_numerical': 0.011095779067349678, 'mean_word_length': 4.210676647296365, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Z2 topological order and first-order quantum phase transitions in systems with combinatorial gauge symmetry

We study a generalization of the two-dimensional transverse-field Ising model, combining both ferromagnetic and antiferromagnetic two-body interactions, that hosts exact global and local Z2 gauge symmetries. Using exact diagonalization and stochastic series expansion quantum Monte Carlo methods, we confirm the existence of the topological phase in line with previous theoretical predictions. Our simulation results show that the transition between the confined topological phase and the deconfined paramagnetic phase is of first-order, in contrast to the conventional Z2 lattice gauge model in which the transition maps onto that of the standard Ising model and is continuous. We further generalize the model by replacing the transverse field on the gauge spins with a ferromagnetic XX interaction while keeping the local gauge symmetry intact. We find that the Z2 topological phase remains stable, while the paramagnetic phase is replaced by a ferromagnetic phase. The topological-ferromagnetic quantum phase transition is also of first-order. For both models, we discuss the low-energy spinon and vison excitations of the topological phase and their avoided level crossings associated with the first-order quantum phase transitions.


I. INTRODUCTION
Topological quantum states of matter are of central focus in modern condensed matter physics. One of the main features of strongly-interacting systems with gapped topological order is that they can present degenerate ground states. This degeneracy cannot be lifted by the action of local perturbations, and hence this property makes such systems perfect candidates for building stable (topological) qubits. Several theoretical models have been proposed to realize gapped topologically ordered states. For instance, the toric code [1] and dimer models on non-bipartite lattices [2,3] host quantum spin liquid (QSL) phases that possess Z 2 topological order. In both these examples, the Hamiltonians contain multi-body interactions, making it a challenge to encounter materials realizing these phases or to realize them in artificial structures.
Attempts have been made to construct models with simpler interactions that can host gapped QSLs. A rare successful example is the cluster charging model of bosons on the kagome lattice [4][5][6], which has been shown theoretically and numerically to host a Z 2 quantum spin liquid, in a system with only two-body interactions. However, these two-body interactions are of the XXZ type, which are not easily implementable in, say, programmable quantum devices. Moreover, the Z 2 gauge symmetry in this model is only emerging, i.e., it exists in the effective model derived in perturbation theory, but it is not an exact symmetry of the original Hamiltonian.
Recently, a construction for which the Z 2 gauge symmetry is exact was proposed on a variant of the transverse field Ising model (TFIM), utilizing only simple two-body ferromagnetic and antiferromagnetic ZZ interactions [7]. Monomial (matrix) transformations that correspond to combinations of spin flips and permutations play a central role in the construction, thus dubbed combinatorial gauge symmetry. Because the construction utilizes only ZZ interactions (of both signs) and a transverse field, the model can be easily implemented, for example, with current Noisy Intermediate-Scale Quantum (NISQ) technology using flux-based superconducting qubits, or other types of quantum computer architectures that provide similar interactions on qubits. The model has already been successfully implemented on a D-wave quantum device in a recent experiment [8].
In this paper, we present a quantitative and detailed study of the combinatorial Z 2 gauge model originally proposed in Ref. 7. Two different types of quantum fluctuations are introduced while preserving the gauge symmetry: a transverse field acting on the gauge spins and a XX ferromagnetic interaction between the gauge spins, respectively. In both cases, we observe the existence of a Z 2 topological state separated by a first-order transition from another phase-a paramagnet in the first model and a ferromagnet in the second The structure of the paper is as follows. First, in Sec. II we give a brief introduction to the model with combinatorial Z 2 gauge symmetry. The model realizes the Z 2 gauge symmetry through monomial transformations and effectively realizes the 4-body interaction term as the star term in the classical version of the toric code. We further introduce two types of quantum fluctuations by applying either a transverse field on the gauge spins (model-X) or a XX ferromagnetic interaction between the gauge spins (model-XX), both of which respect the gauge symmetry. In Secs. III and IV, we provide numerical results on both models obtained from quantum Monte-Carlo (QMC) simulations with the Stochastic Series Expansion (SSE) method as well as exact diagonalization (ED). In both cases, we find that the system exhibits a Z 2 topologically ordered phase separated by a first-order transition from either a paramagnetic phase (model-X) or a ferromagnetic phase (model-XX). We summarize our main results and discuss the remaining open questions and future prospects in Sec. V.

II. SPIN MODELS WITH COMBINATORIAL GAUGE SYMMETRY
We consider two models in our studies starting from a baseline model that contains ZZ interactions between the spins with both ferromagnetic and antiferromagnetic couplings arranged in a pattern such that the combinatorial gauge symmetry is realized.
On top of the baseline model, additional kinetic terms are introduced to give the system quantum dynamics. Here we consider two different types of kinetic terms. The first one is a transverse field on the gauge spins in what we call model-X; the second one is a ferromagnetic XX-type (σ x σ x ) coupling between the gauge spins, which defines model-XX. In this section we start by describing the details of the baseline model and how the combinatorial gauge symmetry is realized. We then discuss the models with the two different kinetic terms.

A. Baseline model with combinatorial gauge symmetry
In the baseline model S = 1/2 spins reside on both the sites and the links of a square lattice, as shown in Fig. 1(a). For each star (or vertex) s of the lattice, we place four "matter" spins µ (the orange dots inside the circles representing the lattice sites) and four "gauge" spins σ (green square) on the links of the star. One such star with its total of eight spins is depicted in Fig. 1(b) along with a labeling scheme. The Hamiltonian of the system is written as where H σ kin is the kinetic term involving only the σ x component of the gauge spins (on all links), and at each star s we define a local Hamiltonian on its eight spins Here the sums with j and a indices include the four matter spins σ z j and the four gauge spins µ z a , respectively. Notice that H s contains a transverse field only on the matter spins. The ZZ couplings between the gauge and matter spins have magnitude J and signs controlled by the Hadamard matrix W , leave the spin commutation relations invariant if L and R are monomial matrices, i.e., generalized permutation matrices with a ±1 entry in each line or column in the case of the group Z 2 . The transformations correspond to combinations of rotations by 0 (+1 entry) or π (-1 entry) around the x-axis, followed by a permutation of the indices. A local gauge symmetry is generated by flipping gauge spins on closed loops around the elementary plaquettes, together with accompanying transformations on matter spins. Flipping the gauge spins around the loop corresponds to choosing R matrices for each s traversed, with an even number of −1 entries associated with the links visited. For each such R, there is a corresponding monomial matrix L = W R W −1 [7]. These pairs of monomial R and L matrices are such that W = L −1 W R, and thus the transformation in Eq. (4) leaves the ZZ part of the Hamiltonian invariant. Moreover, since in H s the transverse field on the µ x is the same on all a ∈ s, the permutation action of the monomial L also leaves these terms unchanged. Hence, H s is invariant under the monomial transformation with L and R. Finally, since R is diagonal and only σ x enters in H σ kin , this kinetic term is also invariant. For each such loop, we have a Z 2 symmetry of the Hamiltonian in Eq. (1).
This local Z 2 gauge symmetry is exact for any value of the parameters in the Hamiltonian Eq. (1). We can obtain further intuition by connecting to the more familiar formulation of the Z 2 gauge theory [9,10] in certain limits. Consider the effective Hamiltonian for the H s terms when their energy scales are larger than those in H σ kin ; in this regime, one can diagonalize H s by fixing the σ z i around the star and treating the problem as that of a paramagnet for the matter spins µ a . The result is an effective Hamiltonian for the lowest states that take the form of a four-spin interaction among the gauge spins: where the parameters γ and λ are given by [7] Notice that the effective H eff s is, up to a constant shift, the same as the star term A z s = i∈s σ z i that appears in the toric code [1] and the Z 2 lattice gauge model [9,10]. We depict in Fig. 2(a) the star term in our model, juxtaposed to the star term represented in the toric and Z 2 gauge models in Fig. 2(b) as the product of four spins on the red cross. The manifold of other states in our model, those beyond the effective term, are separated by a scale Γ m . Thus, in the limit Γ m → ∞, higher energy sectors are projected out, and the system Hamiltonian asymptotically becomes the exact star term of the toric code.
In the conventional Z 2 gauge transformation, the gauge operator G toric p = i∈p σ x i is built as a product of σ x i spins along the smallest loops, shown as the blue square in Fig. 2(b)). In a system with Hamiltonian as in Eq. (1), which possesses combinatorial Z 2 gauge symmetry, the exact local gauge transformation on a plaquette p includes additional transformations corresponding to the action of an operator L µ s on the matter spins of star s as which implements the flips and permutations associated to the monomial matrix L. The plaquette term generating the combinatorial Z 2 gauge symmetry is then defined as In a system with linear size L (total spins N = L × L × 6) and periodic boundary condition, one can find M = 2 + (L 2 − 1) independent gauge operators G p that commute with the Hamiltonian Eq. (1). Within the M operators, two of them, G x and G y , are defined along non-contractible loops in the two spatial directions, and their quantum numbers uniquely characterize the topological ground state degeneracy. Therefore in the basis of these G p operators G p |q = q |q and the Hamiltonian can be block-diagonalized into 2 M blocks associated with unique quantum number sets q = (G x = ±1, G y = ±1; ±1, ±1, . . . ) (See Appendix. A for details), where we list the quantum numbers of the non-contractible loops G x and G y first, and the remaining quantum numbers are associate to the other M − 2 independent local gauge operators. The discussion thus far is rather general, and, in particular, the combinatorial gauge symmetry is exact, provided that the kinetic term H σ kin involves only the σ x component. Below we shall discuss two different choices of H σ kin .
B. Model with transverse field on the gauge spins-model-X A simple choice of kinetic term is to apply a transverse field on the gauge spins, or, equivalently, the case with full Hamiltonian In the limit Γ m ≫ J, we can replace the first two terms by the star equivalent Eq. (5a). Therefore, Hamiltonian Eq. (9), which obeys the exact local combinatorial Z 2 gauge symmetry, has the usual Z 2 lattice gauge model as its low energy description: with g = Γ g and λ given by Eq. (5c). The Z 2 lattice gauge model has been well studied and shown to host Z 2 topological order for g < g c , and undergo a deconfinementconfinement transition at g = g c that leads to a paramagnetic phase for g > g c . The critical point is g c ≈ 0.3285λ from the exact mapping to the dual 2D Ising model [11][12][13][14]. We expect similar phases to appear in the full Hamiltonian Eq (9) as the uniform transverse field on the gauge spins is varied. The limit Γ m ≫ J therefore allows us to estimate the phase boundary of our model Eq (9) perturbatively by relating the parameters Γ m , Γ g , and J to the couplings g and λ of the Z 2 lattice gauge model using Eq. (5a)) and Eq. (10) as Setting g = g c ≈ 0.3285λ results in the leading-order phase boundary in the (Γ g , Γ m ) plane. By comparing with unbiased numerical results, we will later show that this Γ m → ∞ form provides a good approximation to the phase boundary even down to values as small as Γ m ≈ 2 [ Fig. 4(a)]. Note that the mapping to the simpler lattice gauge model is only perturbatively exact in the limit of Γ m → ∞. In the region with a small Γ m , we cannot in general rule out other phases appearing away from the perturbative regime. However, we do not see any obvious reasons to expect phases beyond those of the effective Z 2 lattice gauge model. We present in Sec. III our numerical results that support the absence of other phases in our model and reveal the order of the phase transitions between the topological and paramagnetic phases. Another simple choice of kinetic term is to add twospin XX interactions between nearest-neighbor gauge spins, so that the full Hamiltonian is Notice that H σ kin satisfies the general conditions presented above to retain the combinatorial gauge symmetry. (H σ kin is written in terms of σ x only, and thus commute with the local operators G p .) To gain intuition about the model with Hamiltonian Eq. (13), we again consider the limit Γ m ≫ J in which the first two terms can be replaced by their equivalent star term Eq. (5a). Defining the four-spin star operator A z s ≡ i∈s σ z i , the effective Hamiltonian, without the kinetic term, reads H eff 0 = −λ s A z s . The ground states of H eff have A z s = 1 (parity +1) on all stars or vertices. Similarly to the dual mapping of the conventional Z 2 lattice gauge model, we can introduce a conjugate star operator A x s that flips the eigenvalue of A z s , in terms of which we write the gauge spin between two nearest-neighbor stars s and s ′ as σ x i = A x s A x s ′ . The XX kinetic term Eq. (12) perturbatively generates four-spin plaquette interactions (i.e., products of four σ x around the small loop around a plaquette). To this end, the bond operators in the kinetic term can be arranged, e.g., in the way illustrated in Fig. 3. Starting from the ground state of H eff 0 , where all vertices have parity A z s = +1, acting with a term of the form σ x k σ x l on two sites k and l within a vertex s generates a pair of defects on either the A or the B sublattice, depending on the bond chosen. Adding another bond operator within the same plaquette, parallel to the first bond, generates a plaquette term at second order in J x /λ. The effective perturbative Hamiltonian in terms of the A z,x operators then takes the form where ss ′ indicates nearest-neighbor stars within the same sublattice, and λ is the aforementioned effective strength of the star-term, Eq. (5c). In the limit Γ m ≫ J, our model effectively reduces to two independent TFIMs; one on each sublattice. The particular arrangement of the XX terms imposes an additional even-odd conservation law in the system. The parity of negative A z s vertices within a sublattice is conserved, which can be easily seen from Fig. 3(a), where the XX operators can only create defects in pairs within each one of the sublattices.
In the limit J x → ∞, since all gauge spins interact ferromagnetically, an x-direction ferromagnetic phase arises; thus we expect a quantum phase transition between the topological phase and a ferromagnetic phase, replacing the topological-paramagnetic transition of the model-X discussed in Sec. II B. Similarly to the model with transverse field on the gauge spins, where the phase boundary between the topological and paramagnetic phases is given perturbatively by Eq. (11) through the mapping to the TFIM, here this mapping gives the following relation between the field strength in the TFIM and the field λ in Eq. (5c); Thus, the topological-ferromagnetic phase boundary can be obtained to leading order by setting the ratio above to the critical point 3.04497(18) of the 2D quantum Ising model [12,13]. This boundary also describes numerical results for surprisingly small values of the matter field, down to γ m ≈ 1 [ Fig. 11(a)].

III. ANALYSIS OF MODEL-X
In this section, we present our numerical studies of the model-X introduced in Sec. II B. Our results support the theoretical conjecture that the model has a topological and a paramagnetic phase with no other phases. The nature of the quantum phase transitions between these two states is revealed. In the following, we fix J = 1 and impose periodic boundary conditions in all our numerical simulations.

A. Fidelity susceptibility
To confirm that the model does have phases predicted by the effective Z 2 gauge theory, we start by identifying signatures of the phase transition. We first consider the fidelity susceptibility which can probe the existence of a phase transition without requiring knowledge of any order parameter [15]. The fidelity susceptibility is defined as the second derivative of the logarithmic fidelity with respect to a generic tuning parameter x where | is the infinitesimal fidelity in the direction defined by x. In our model with transverse fields, two types of fidelity susceptibilities can be formally defined with variations along the two different transverse fields; x = Γ m or x = Γ g ; Across a phase transition, the fidelity susceptibility should develop a maximum that diverges in the thermodynamic limit. We first calculate the fidelity susceptibilities using exact diagonalization with the Lanczos method [16,17] on a small system with 2×2 unit cells, i.e., N = 24 spins in total. Due to the rapid growth of the Hilbert space, this is currently the largest accessible system size with our computational resources. Fig. 4(b) shows χ m F as a function of Γ m for several different values of Γ g . A single peak is present in all cases, which implies the possibility of a phase transition in the thermodynamic limit. Furthermore, we have not observed any cases of multiple maxima in any of our calculations, suggesting only two different phases. Similar behaviors are also observed for χ g F as shown in Fig. 4(c). From the location of the maximum of χ m F shown in Fig. 4(a), we find that the data fall close to the perturbative (large-Γ m ) topological-paramagnetic phase boundary even though the value of Γ m is not extremely large (and Γ g not extremely small).
It may seem surprising that the phase boundary is given accurately by a system with only four unit cells. To confirm that the observed maximum of the fidelity susceptibility grows with the system size and truly indicates a phase transition, we next turn to QMC simulations to reach larger system sizes. We use the SSE QMC method [18,19], for which a convenient way to compute the fidelity susceptibility was devised recently [20]. In our simulation, we set the inverse temperature as β = 4L. Because of the small vison gaps in the topological phase, this scaling of β does not allow us to reach the finite-size ground state deep inside the topological phase. However, with the T → 0 limit approached with β ∝ L we can still address the nature of the quantum phase transition from the gapped paramagnetic phase. In our model, the Ising interactions are highly frustrated, and to mitigate the associated effects of slow dynamic of the QMC updates in the topological phase and at the phase transition, we have implemented quantum replica exchange [21][22][23]. Simulations are thus carried out in parallel for a large number of replicas with different values of Γ m on both sides of the transition, with swap attempts carried out for neighboring values of the parameter after several conventional SSE updates. Even with replica exchange, it is still difficult to equilibrate systems for large L, and we have limited the present study to L ≤ 8. As we will see, these moderate system sizes are already sufficient for drawing definite conclusions. Fig. 5(a) shows the results of the fidelity susceptibility χ m F at Γ g = 0.2 as a function of Γ m . We indeed find that the peak identified in the ED calculations diverges upon increasing the system size, providing solid evidence of a phase transition. The other fidelity susceptibility χ g F shows a similar behavior, as shown in Fig. 5(b).
In order to understand the nature of the phase transition, we perform finite-size scaling of the maximum value of the fidelity susceptibility to extract the associated critical exponent. Based on the similarity of the model to the Z 2 lattice gauge model, in which the transition is in the (2+1)D Ising universality class, one might naively expect a continuous transition at which the maximum should scale with system size as [16] χ F /L d ∼ L 2/ν−d with d = 2 the spatial dimensionality and ν ≈ 0.63. However, we do not observe a scaling of the above form. Instead, we analyze the data using a generic scaling form To further take into account finite-size corrections, we consider a size dependent exponent b L extracted from two system sizes, L and L + 2; graphically this exponent corresponds to the slope of the line drawn between two data points on a log-log scale as shown in the insets of Fig. 5.
In the case of χ m F in Fig. 5(a) we find b 4 ≈ 3.0 (i.e., the line drawn between data points for L = 4 and L = 6) and b 6 ≈ 3.5 (from L = 6, 8). In the case of χ g F we find b 4 ≈ 2.7 and b 6 ≈ 3.2. These exponents are significantly larger than the expected value 2/ν − d ≈ 2/0.63 − 2 ≈ 1.175 of the (2+1)D Ising universality class, and for both susceptibilities the deviation becomes larger for the larger system sizes. It therefore appears more likely that the transition is first-order. Generally, at classical first-order transitions the same scaling forms hold as for continuous transitions, but with the exponent 1/ν replaced by the dimensionality d [24,25]. In a quantum system, the replacement should be 1/ν → d+z, where the appropriate value of the dynamical exponent z reflects the nature of the low-energy excitations in the two coexisting phases [26,27]. Our results in Fig. 5 suggest a first-order behavior with z = 1, in which case We show this type of divergence for reference with the dotted lines in the insets of Figs. 5(a) and 5(b); this asymptotic behavior seems very plausible based on the available data.

B. Topological order
Next, we turn to the properties of the underlying phases. Based on the mapping to the Z 2 lattice gauge model, we expect the phase with small Γ m to be a Z 2 topological quantum spin liquid. Note that Elitzur's theorem forbids any spontaneous symmetry breaking of local gauge symmetries; thus one cannot define any lo- cal order parameter to characterize such topological order [10,28,29]. To detect the topological order, we investigate the global, non-contractible Wilson loop operator, defined as the product of gauge spin σ z operators along a non-contractible loop in the α-direction for the set of sites {i} α,n belonging to the n-th row or column. For a Z 2 spin liquid, the quantum numbers p x,n = ±1, p y,n = ±1 characterize the four degenerate (in the thermodynamic limit) topological ground states regardless of which row or column n is chosen. We can take advantage of this property to define a correlation func- tion detecting the topological order using the product of two parallel non-contractible loops on rows or columns labeled by m and n: which we also illustrate in the inset of Fig. 6. Instead of investigating this correlation as a function of the distance r mn between the two loops, we here take the longest distance for a given lattice size, r mn = L/2, and analyze the L dependence of C p (L/2) defined as a summation over all translations (to reduce the statistical fluctuations) of the two Wilson loops oriented in the α ∈ {x, y} lattice direction: which can be averaged over the two directions. In the topologically ordered phase we expect C p (L/2) = 0 when L → ∞, while in the paramagnetic phase C p (L/2) → 0. Note that this quantity has been used in a previous study of topological order in classical Ising gauge models at zero and non-zero temperatures [30]. In Fig. 6, we show SSE results at Γ g = 0.2 as a function of Γ m . We see that C p (L/2) indeed vanishes with increasing L for large Γ m , while it converges to a finite value for Γ m in a range consistent with the transition point found above for the same value of Γ g . Below we will discuss the size-extrapolated phase boundary. We stress here that the Wilson loop order parameter does not detect any phases with only local order parameters (See Appendix. D for a study of the ferromagnetic state as an example) and our results therefore demonstrate conclusively a Z 2 topological phase of finite extent as the field Γ m is turned on.
Having established a good topological order parameter, we further provide evidence of a first-order phase transition by a finite-size scaling analysis of the corresponding Binder ratio. To this end, we define the topological order parameter on the entire system as the sum of Wilson loops with P 2 = P 2 x + P 2 y where with α ∈ {x, y}, and the Binder ratio In a perfect Z 2 ordered topological phase, P x = ±1 and P y = ±1, forming a Z 2 × Z 2 symmetric order parameter distribution, while P x = P y = 0 in the paramagnetic phase. With increasing system sizes, the Binder ratio is expected to form a step function at the transition point in the case of a continuous transition, while the distribution of the order parameter in the coexistence state at a first-order transition typically is also associated with a divergent peak adjacent to the step [31,32]. In Fig. 7(a) the Binder ratio indeed evolves into a step function with a side peak, though the latter is only seen clearly for the largest system sizes, L = 8, and for L = 6 there is a very weak maximum as well. Looking at the derivative of B, in Fig. 7(b) we observe a divergent positive main peak, and for L = 8 a prominent negative peak reflects the presence of the first-order side peak in Fig. 7(a). Thus, we have strong evidence of phase coexistence at a firstorder quantum phase transition caused by an avoided level crossing. In Appendix E, we further show the full distribution of the two-component Wilson loop order parameter ρ(P x , P y ), which clearly shows a phase coexistence characteristic of a first-order transition. The maximum of the Binder ratio derivative diverges with the system size as the step function develops. The derivatives can be evaluated directly in the SSE simulations, using an estimator derived in Appendix B. With the replacement 1/ν → d + z in the scaling form dB/dΓ m ∼ L 1/ν and expecting z = 1 here, the peaks should diverge as L 3 . Indeed, in Fig. 7(b) the peak values for the three largest system sizes can be fitted to a power-law L b with b ≈ 3.0, thus supporting a first-order transition.
So far, we have discussed the divergence properties of the peaks in the Binder ratio and the fidelity susceptibility. We also need to extrapolate the peak locations in order to obtain the transition point in the thermodynamic limit. Fig. 8 shows the dependence of the peak locations on 1/L along with extrapolations assuming power-law corrections. All quantities show mutually consistent behaviors for the largest system sizes, but χ m F has much larger scaling corrections than the other quantities. Extrapolations with error analysis give the critical value of the matter field Γ c m ≈ 2.27 for the gauge field Γ g = 0.2 considered here. In the phase diagram in Fig. 4(a) we have marked this transition point with a circle, and we also show the result obtained using the same methods for Γ g = 1. These QMC points are very close to the boundary estimated from the ED results for a very small system with L = 2.
Here we should note that the ED results are calculated exactly at T = 0, while there are still some temperature effects left in the QMC results obtained with our choice of temperature scaling, T = (4L) −1 . In the case of L = 2, the QMC results for Γ c m (L) are actually quite far from the T = 0 ED results because of the temperature effects. However, since T → 0 as L increases, the L → ∞ extrapolated QMC results are not affected by finite temperature (beside unimportant constant factors in the peaks of the physical quantities). In this regard, it can also be noted that effects of inappropriate temperature scaling with L could potentially ruin a quantum phase transition that does not extend to T > 0 (as is the case with topological order in two spatial dimensions), while there is no reason to expect a transition detected when T ∝ 1/L to vanish if T approaches zero more rapidly.

C. Energy derivative
In this section, we show another signature of the firstorder phase transition from the ground state energy density. In practice, the T → 0 internal energy is obtained with QMC calculations at T = (4L) −1 . The type of first-order transition indicated by the results above, where finite-size scaling with the exponent replacement 1/ν → d + z holds, should be associated with an avoided level crossing. Thus we expect a change in the slope of the energy at the transition with increasing L. In Fig. 9(a) we plot the energy per site as a function of Γ m for the same gauge-field strength as considered above, Γ g = 0.2. At first sight, the data exhibit a smooth behavior without any visible kink. However, by taking the numerical derivative of the energy, as shown in Fig. 9(b), we find a clear signature of non-analytic behavior, such that the derivative becomes discontinuous at the transition in the thermodynamic limit.
Along with the other results, this demonstration of a discontinuous energy derivative provides definite proof of a first-order quantum phase transition between the topological and paramagnetic states.

D. Level spectroscopy
Having used QMC simulations to establish the existence of an extended Z 2 topological phase and its quantum phase transition into the paramagnetic phase, we now again turn to Lanczos ED calculations in order to investigate the energy level spectrum of the system. We use the combinatorial Z 2 symmetry to block-diagonalize the Hamiltonian into M = 2+L 2 −1 blocks in the basis of Z 2 gauge generators, as discussed in detail in Appendix A. The blocks are categorized by a set q = (±, ±, {±}) of quantum numbers, where the first two elements correspond to the two non-contractible loops with associated quantum numbers G x and G y and {±} denotes the set of L 2 − 1 local quantum numbers G i . In the thermodynamic limit, the topological ground state should be four-fold degenerate, corresponding to the lowest energy states from sectors with G x = ±, G y = ± and G i = + for all other i ∈ M − 2 local operators. The transverse field does not commute with the Hamiltonian; thus, there are always finite-size gaps between the four topological states in a finite system. Our ED calculations here are again restricted to L = 2 (and we present some QMC results also for L = 4), but even for this very small system many of the salient signatures of spinon and vison excitations FIG. 10. ED level spectrum at Γg = 0.2 for a system of size L = 2 (N = 2 × 2 × 6 spins). In (a) the two lowest energy gaps relative to the ground state (marked Gnd in the legends) are graphed versus the matter field strength for each of the 32 blocks with quantum numbers q = (±, ±, ±, ±, ±). Many blocks are degenerate because of lattice symmetries; thus the number of different curves is much less than 64. The four states that become degenerate topological ground states in the Z2 phase are marked by blue crosses (the finite-size ground state), red stars, and green triangles (two degenerate sectors); these states all have the local quantum numbers Gi = +1. The blue dots represent all other lowest block levels; these are states with visons (two or a larger even number) of the topological phase (marked by v in the inset illustration). The orange curves represent the second-lowest states in each block; they correspond to the spinon excitations (particles indicated in the inset by s) of the topological phase and they all become degenerate for Γm → 0. In (b), the lowest block states are graphed on a magnified scale. The eight-fold degenerate level that is the lowest excitation in the paramagnetic phase was calculated with both ED (yellow dots) and extracted from imaginary-time correlations from QMC simulations (violet circles), to demonstrate the correctness of the latter for L = 2. QMC results for both L = 2 and L = 4 (calculated at inverse temperature β = 24L) are shown in the inset.
can be observed, as well as signatures of the first-order quantum phase transition.
In Fig. 10(a), we graph low-energy gaps ∆ E relative to the ground state versus the matter field Γ m at fixed Γ g = 0.2. For each of the 32 topological symmetry blocks pf the L = 2 system, the two smallest gaps are shown, but because of degeneracies due to lattice symmetries there are only 11 distinct curves. The unique finite-size ground state has q = (+, +, +, +, +), i.e., G x = +, G y = + and G i = + for i = 1, 2, 3. The four levels that become degenerate in the thermodynamic limit in the topological phase are highlighted with different symbols. As for the remaining low-energy levels, we note that in the Z 2 topological phase two types of quasi-particle excitations should be expected; spinons (s) and visons (v), which are created in pairs by acting on the ground state with σ x and σ z respectively on the gauge spins, as indicated in the insets of Fig. 10(a).
In the Γ m → 0 limit, the spinon excitations are gapped with ∆ E ≈ 4J (the gap value is 4J if Γ g = 0), as seen clearly in Fig. 10(a), where these levels are shown with orange symbols. The vison gap opens when increasing Γ m , as can be seen from the fact that the effective model takes the form to lowest order in perturbation theory. Here G p is the local gauge generator that appears at 12th order, where the coupling is We can identify the vison excited states simply by considering the quantum number blocks that couple to the ground state through the on-site σ z operators (which do not commute with the local gauge operators). In Fig. 10(b), we plot ED results for the same parameters as in Fig. 10(a), but with a change in scale to focus on the vison states. These states are gapped for all Γ m > 0, but the gaps are much smaller in the topological phase than in the paramagnetic phase. The lowest vison state, which contains two visons, is eight-fold degenerate on the small system considered here. The other levels in Fig. 10(b) correspond to states with (an even number) more than two visons. We can extract the lowest vison gap from QMC simulation by analyzing the imaginary-time autocorrelation function of σ z , defined as where N g = L × L × 2 is the total number of gauge spins in the system. The gap can be extracted by fitting an exponential function to G(τ ) for large τ . Here it is also important that the temperature is sufficiently low, so that the asymptotic form of G(τ ) is dominated by the lowest gap; See Appendix. C for further technical details on these calculations. In Fig. 10(b) we compare the lowest gap extracted from the QMC data for the L = 2 system at inverse temperature β = 48 with the ED result. We observe good agreement between the two calculations. Note that the eight-fold degenerate levels with two vison excitations undergo a true (not avoided) level crossing with a state with G x = G y = + and all local quantum numbers G i = −. This level-crossing is a finite-size effect, and we do not expect such behavior to persist for larger system sizes. The G i = − state contains four visons, i.e., it can be reached from the ground state with two different σ z i operations. It therefore does not contaminate the correlation function G(τ ) corresponding to the two-vison level of interest. For this small system the four-vison state falls under the lowest 2-vison state below Γ m ≈ 3, i.e., close to the phase transition into the topological state.
Note again that the quantum numbers G i are conserved (i.e. commute with the Hamiltonian) in both the topological phase and the paramagnetic phase of the model. However, visons with G i = − are deconfined only in the topological phase. In the paramagnetic phase, the lowest energy vison excitations must form bound states residing on two adjacent plaquettes, while states with more visons and larger separations between the visons have larger energy costs, as shown in Fig. 10(b).
In the inset of Fig. 10(b), we show the vison gap based on QMC calculations for both L = 2 and L = 4 at Γ m = 0.2, using inverse temperature β = 24L. While the L = 2 gap exhibits only a rather smooth variation with Γ m , at L = 4 a sharp feature has developed close to the phase transition. The sharp behavior of the gap here is consistent with the scenario of a first-order transition through an avoided ground state level crossing, and this mechanism should be associated also with avoided level crossings of the low-lying excitations.
Although the model possess Z 2 topological order and can be directly implemented on existing quantum devices [8], it will be difficult to reach the true ground state, or even a thermal state with a low density of visons, due to the very small vison gap. These difficulties are clear from the effective model obtained perturbatively, Eq. (23) with the 12th-order effective coupling in Eq. (24). Nonetheless, there may still be signatures of the mutual statistics of the spinons and visons that could be observed in the regime where temperature is larger than the vison gap but still much smaller than the spinon gap, as discussed in Ref. [33]. In this regime the visons randomly appear within plaquettes because their energy of formation is much smaller than the temperature. In the presence of kinetic terms (such as a transverse field), the spinons acquire dynamics at a scale much faster than that of the visons, so effectively they quantum diffuse in a background or randomly placed visons. Because of the mutual statistical phase of π between the two types of particles, the random visons serve as sources of π fluxes, which lead to quantum interference corrections to the diffusion of the spinons.

IV. ANALYSIS OF MODEL-XX
Here we present numerical results for the model-XX introduced in Sec. II C, organized in the same way as the numerical studies of model-X in Sec. III. Lanczos ED calculations were carried out for L = 2 systems. QMC results for larger systems up to L = 6 were again obtained using the SSE method supplemented with quantum replica exchange. We fixed Γ m = 1.0 and simulated several replicas at different values of J x across the two phases. The values are chosen such that the acceptance rate of swapping neighboring replicas is in the range 0.4 ∼ 0.6. Our findings and arguments are very similar to those for model-X discussed in Sec. III, with the exception of issues pertaining to the ferromagnetic phase, and we therefore keep the discussion brief in this section.

A. Fidelity susceptibility
As in the previous study of model-X, we first discuss Lanczos ED results for the fidelity susceptibility χ x F defined as in Eq. (17b) with the substitution Γ g → J x . In Fig. 11(b), we show χ x F versus J x for several values of Γ m . In all cases, we observe a peak indicative of a phase transition. The locations of the maxima are shown along with the perturbative phase boundary in Fig. 11(a).
In Fig. 12(a) we plot SSE results for larger systems. As expected we find a maximum that diverges with increasing system size, showing a true phase boundary and only two phases. Analyzing the peak height using the effective exponent b L defined with system sizes L and L + 2, we find b 2 ≈ 4.8 and b 4 ≈ 3.9. These exponents are again significantly larger than the 2/ν − d ≈ 2/0.63 − 2 ≈ 1.175 for the (2+1)D Ising universality class, but close to b = 2(d + z) − d = 4 for a first-order transition when z = 1. We note one difference with respect to the previous model, as seen in Fig. 5, in that case the exponent b L increases with L, while in the present case it decreases. One may then question whether the value b = 4 is really obtained in the limit L → ∞. Nevertheless, all our complementary results to be presented below also lend support to a first-order transition.

B. Topological order
To detect the topological order, we again consider the correlation function C p (L/2) between two parallel noncontractible Wilson loops, defined previously in Eq. (20). Fig. 12(b) shows results at Γ m = 1.0. Here a discontinuity reflecting the first-order transition develops more clearly as compared to the results for model-X in Fig. 6, thus suggesting a more strongly first-order transition in this case. Note, however, that the parameter values chosen for the two models in these figures, Γ g = 0.2 and Γ m = 1, are not directly comparable. In both cases, the strength of the discontinuity will vary with the model parameters.
As in Sec. III, we use the Wilson loop order parameters P x = a p x,a and P y = a p y,a to confirm the extent of the topological phase. We use the Binder ratio B as defined in Eq. (22) with both components taken into account. The results, shown in Fig. 13(a) exhibit developing step functions with associated peaks indicative of a first-order transition. The derivatives exhibit the expected divergent peaks. Because of the limited system sizes, we refrain from analyzing the peaks further. We have used the peaks to extrapolate the transition point to infinite size and show the result with the green circle in the phase diagram in Fig. 11 at Γ m = 1. As in the model-X with transverse field on the gauge spins, we find only a small difference between the QMC result and the L = 2 ED result in this case.

C. Energy derivative
We present further evidence of a first-order phase transition from the energy density. As shown in Fig. 14(a), in this case we observe a clear kink behavior for the larger system sizes, L ≥ 4, and the derivative in Fig. 14(b) accordingly shows a strong discontinuity developing.

D. Level spectroscopy
We have again used Lanczos ED to find low-lying states for each block of quantum numbers characterizing the combinatorial Z 2 symmetries in the L = 2 system. In Fig. 15, we present the two smallest gaps versus the XX coupling J x at Γ m = 1.0. The lowest states in the sectors G x = ±, G y = ± and G i = + again are those that will eventually become degenerate as L → ∞ in the topological phase, and these states are highlighted with different x + P 2 y . The divergent peak next to the step indicates phase coexistence at a first-order transition. Accordingly, in (b) the derivative of the Binder ratio shows divergent positive and negative peaks.
symbols. The lowest energy excitations in the topological phase, states with visons, form levels very similar to what we saw in the model-X. However, the spinon spectrum looks very different. Due to the XX ferromagnetic interaction, spinons created in pairs within one of the sublattices has lower energy comparing to the one created in neighboring pairs (created by a single σ x operation) as illustrated in Fig. 15(a).
At J x → 0, only spinon excitations exist, with a gap size of order O(J), as can be seen in Fig. 15(a) where these levels are marked in orange. The vison gap opens with increasing J x , as the effective model from the lowest order in perturbation takes the form Here G p is the local gauge generator (plaquette term) that appears at 10th order, with g ∼ Γ 8 m J 2 x /J 9 , which should be compared to the 12th order perturbative Hamiltonian in the case of model-X.
In the case of the XX interaction used here, there is an additional gauge spin inversion symmetry in x-basis that is not present in the model-X. Define the inversion operator V = i σ z i as the product of all σ z gauge-spins. This operator clearly commute with Hamiltonian and its quantum numbers v = ± correspond to symmetric or antisymmetric states. We find that all the lowest energy levels of the 32 gauge blocks (blue) are symmetric and the second state (orange) is always anti-symmetric except for the highest energy level shown in Fig. 15(a), which exhibits an actual level crossing (see appendix F for further details).
Among all the spinon excitations, the lowest one belongs to the same sector as the ground state, with q = + + + + +. This excitation, which is marked with a dashed line in Fig. 15(a), becomes degenerate with the ground state for large J x , reflecting the ferromagnetic Ising order with spontaneously broken Z 2 symmetry in the thermodynamic limit.
For the vison excitations, since σ z operators do not commute with local gauge operators, we can identify the vison excited states simply by considering the quantum number blocks that couple to the ground state through the on-site σ z operators. In Fig. 15(b), we plot ED results for the same parameters as in Fig. 15(a), but with a change in scale to focus on the vison states. These states are gapped for all J x > 0, but the gaps are much smaller in the topological phase than in the FM phase. The lowest vison state, which is marked by yellow symbols in Fig. 15(b), contains two m particles, and is eight-fold degenerate on the small system considered here. The other levels marked by blue in Fig. 15(a) correspond to states FIG. 15. Level spectrum relative to the ground state for the L = 2 system with XX interactions, calculated with Lanczos ED. The organization of panels (a) and (b) is as in the corresponding Fig. 10 for the model-X. We refer to the same for further explanation of the visualization of the spectrum. The lowest spinon excitation, shown as orange circles with a dashed line, is the first excited state from the same block as the ground state (which has q = + + + + +). This state is antisymmetric with respect to spin inversion and forms the two-fold degenerate multiplet together with the corresponding symmetric state in the ferromagnetic phase. with more than two (an even number of) visons.

V. CONCLUSIONS AND DISCUSSION
We have presented a numerical study of spin models with only one-and two-spin interactions that realize a combinatorial Z 2 gauge symmetry. We considered two models that only differ by the kinetic terms given to the gauge spins: model-X (containing a transverse field) and model-XX (containing XX interactions). We found conclusive evidence for an extended Z 2 topological quantum spin liquid phase in both models.
In the case of model-X, we identified two phases; a topological phase and a paramagnetic phase. We demonstrated a first-order quantum phase transition between these phases, in contrast to the well known continuous transition of the conventional Z 2 lattice gauge model. In model-XX we identified a topological phase and a competing ferromagnetic state. Our data also support a firstorder transition between these two phases in model-XX.
Perturbatively, the XX interaction of model-XX generates a plaquette operator G p at a lower order in perturbation theory as compared to the transverse field of model-X, and therefore the size of the vison gap increases, as we also observe.
The presence of the first-order transition between the topological and the competing state, in both models, raises the following interesting question: As we have discussed in the paper, in the limit of a large transverse field Γ m on the matter spins, the models map to the usual Ising gauge model, which has a continuous transition. An important question is then whether the continuous transition persists for some finite range of values of Γ m , or whether it turns first-order immediately. This question can in principle be answered by considering the corrections to the usual Z 2 gauge model in Eq. (10), which will appear when carrying out a perturbative expansion to higher order in Γ −1 m . The question is then whether these corrections are renormalization-group relevant or irrelevant at the critical point. While we have not carried out this expansion and duality mapping, it appears likely that the additional interactions generated in the Ising model will involve products of more than two spins, and most likely these interactions will be irrelevant at the Ising critical point. Thus, we suspect that there will be indeed a tricritical point separating continuous Ising transitions and first-order transitions for large values of Γ m in Figs. 4 and 11. We leave tests of this hypothesis open for future work. In analogy to the standard Z 2 lattice gauge theory, in our model the 4-body interaction term A i is effectively generated with the term (See Ref. [7] for further details) and the monomial transformation leads to the modification of the local gauge generator with the additional L µ s representing the action of monomial matrices L acting on matter spins. All the formulations that characterize the symmetry and quantum numbers in the standard Z 2 lattice gauge theory can also be applied to our model with combinatorial Z 2 gauge symmetry.
There are in total M = 2 + (L x L y − 1) such independent operators in a L x × L y system with periodic boundary conditions. Within the M operators, two of them are defined as a product of σ x , G i ≡ (i) j σ j x along non-contractible loops where their quantum numbers uniquely characterize the 4-fold topological degeneracies of the ground state in the thermodynamic limit. As shown in Figure. 16(a), G y is defined along a noncontractible loop in the y-direction, and G x is defined in the x-direction. Other operators are local, defined as a product of σ x around an elementary plaquette of 4 spins.
Using the fact that G commutes with the Hamiltonian, we can construct eigenstates of the G operators: G |q = q |q where q = ±1. It is straightforward to see that for each G operator, the eigenstate can be constructed by starting from a classical configuration (which we refer to as the "representative" state [18]) via |q = (1 + qG) |rep . Since the Hamiltonian commutes simultaneously with all M operators, the state should be constructed with a product of 1 + q i G i for all G i . As an example, consider a 2 × 2 system as shown in Fig. 16, we where q indicates the quantum number set q = (q x , q y , q 1 , q 2 , q 3 ) = (±, ±, ±, ±, ±), Z is the normalization factor and subscript i indicates the i-th state within the block. There are in total 2 M symmetry blocks in the system.
and we have used the short-hand notation H αiαj = α i |H|α j and O α0α0 = α 0 |O|α 0 . Further, Z is the partition function, M is the operator string length and n is the number of non-identity operators in the current string. The quantity denoted R [H ai,bi ] stands for the product of local Hamiltonian operators, i H αi , where H αi = H i 0 or H i δ . The derivative of the observable with respect to the tuning parameter δ can be calculated from where n δ is the number of the δ operators in the string. We are interested in the Binder ratio of the Wilson loop order parameter, defined as in Eq. (22). Using the above expressions we obtain ∂B ∂δ = ∂ δ P 4 P 2 2 − 2 P 4 P 2 ∂ δ P 2 P 2 4 = 1 δ P 4 n δ + P 4 n δ P 2 2 − 2 P 4 P 2 n δ P 2 To show that the parallel Wilson loop operator does not detect the long range ferromagnetic phase, we perform simulations of the TFIM on L × L square lattices with periodic boundary conditions at inverse temperature β = 2L. We compute the Wilson loop correlator defined in Eq. (20). As shown in Fig. 18, we draw the horizontal line at C P = 1/2 and extract the corresponding value h * x (L) of the transverse field. In the inset of Fig. 18, we demonstrate finite-size scaling of the value h * x as a function of 1/L with a power-law fit of the form y = a + bL −c . Including data for the four largest system sizes, L ≥ 12, we find the best fit with a = −0.02(1), b = 3.15 (2) and c = 0.472 (7), confirming the expectation Wilson loop correlator vanishes for a conventional FM phase in the thermodynamic limit at any finite transverse field. Note that the four peaks at the corners in Fig. 19(a) and Fig. 19(b) should mathematically be of equal size, but they differ here because of the slow migration of the simulation between these peaks, which are separated by tunneling barriers in the SSE configuration space.