Extended Coulomb liquid of paired hardcore boson model on a pyrochlore lattice

There is a growing interest in the $U(1)$ Coulomb liquid in both quantum materials in pyrochlore ice and cluster Mott insulators and cold atom systems. We explore a paired hardcore boson model on a pyrochlore lattice. This model is equivalent to the XYZ spin model that was proposed for rare-earth pyrochlores with"dipole-octupole"doublets. Since this model has no sign problem for quantum Monte Carlo (QMC) simulations in a large parameter regime, we carry out both analytical and QMC calculations. We find that the $U(1)$ Coulomb liquid is quite stable and spans a rather large portion of the phase diagram with boson pairing. Moreover, we numerically find thermodynamic evidence that the boson pairing could induce a possible $\mathbb{Z}_2$ liquid in the vicinity of the phase boundary between Coulomb liquid and $\mathbb{Z}_2$ symmetry-broken phase. Besides the materials' relevance with quantum spin ice, we point to quantum simulation with cold atoms on optical lattices.

There is a growing interest in the U (1) Coulomb liquid in both quantum materials in pyrochlore ice and cluster Mott insulators and cold atom systems. We explore a paired hardcore boson model on a pyrochlore lattice. This model is equivalent to the XYZ spin model that was proposed for rare-earth pyrochlores with "dipole-octupole" doublets. Since this model has no sign problem for quantum Monte Carlo (QMC) simulations in a large parameter regime, we carry out both analytical and QMC calculations. We find that the U (1) Coulomb liquid is quite stable and spans a rather large portion of the phase diagram with boson pairing. Moreover, we numerically find thermodynamic evidence that the boson pairing could induce a possible Z2 liquid in the vicinity of the phase boundary between Coulomb liquid and Z2 symmetry-broken phase. Besides the materials' relevance with quantum spin ice, we point to quantum simulation with cold atoms on optical lattices.
The search of exotic quantum phases with quantum number fractionalization and emergent gauge structure has been an active subject in modern condensed matter physics. One theoretical route in the field is to start from the exotic phase itself and construct solvable models. These models are often contrived and not quite realistic [1][2][3][4] . One exception is the exactly solvable Kitaev model on the honeycomb lattice 5 whose physical relevance to the iridate materials was later pointed out by G. Jackeli and G. Khaliullin 6 . The opposite route is to start from the realistic physical systems and build up relevant models from the physical degrees of freedom. Both routes have been quite fruitful. The latter route faces several major obstacles. Firstly, constructing a relevant physical model itself is not often straight-forward. Secondly, these strongly interacting models often cannot be solved in a controlled manner. Occasionally, certain realistic models, such as the square lattice Heisenberg model for the cuprates, may be solved but yield a bit mundane and known results, and are thus of limited theoretical value for our understanding of strongly correlated quantum matters. Therefore, a physically relevant model, that can be solved in a controlled manner and at the same time gives non-trivial quantum phases, is highly valuable in the study of strongly correlated quantum matters.
The XYZ spin model, that was derived from the microscopics of dipole-octupole doublets on the pyrochlore lattice and on the triangular lattice by one of us and collaborators in Refs. 7 and 8, is a rare example that overcomes the major obstacles of the second route. It was suggested that this model on the pyrochlore lattice could stabilize a U (1) Coulombic liquid 9,10 and may stabilize a Z 2 spin liquid in parts of its phase diagram 7 . The U (1) Coulombic liquid is an exotic quantum state that is described by compact quantum electrodynamics with emergent quasiparticles 9,10 and has found relevance in pyrochlore quantum ice materials 7,11-22 and cluster Mott insulators [23][24][25][26][27] . Besides the non-trivial ground states, it was also pointed out 7 that our model does not have a sign problem for quantum Monte Carlo (QMC) simulation in a large parameter regime, and in fact, it is the case for any lattice 8 . An extension of this model to the kagomé lattice by dimensional reduction from the pyrochlores with magnetic fields was later pursued numerically 28 . Our model was first proposed for various Nd-based pyrochlore materials 7,29-37 , and was recently suggested for a Ce-based pyrochlore spin liquid candidate Ce 2 Sn 2 O 7 38,39 . Thus, the XYZ model becomes a rare model that describes real physical  systems, supports non-trivial quantum phases, and can be solved in a controlled manner in a large parameter regime. Inspired by these compelling properties of the XYZ model 7,8,38 , we carry out both theoretical analysis and numerical calculation to establish the phase diagram of this model on the pyrochlore lattice. We show that the U (1) Coulomb liquid covers a rather large portion of the phase diagram. In addition, the physical boson pairing may render new fates to the emergent spinon-gauge coupling in the U (1) Coulomb liquid 14,15 . We find the thermodynamic evidence for the possible existence of a Z 2 liquid state out of the U (1) Coulomb liquid via an internal Anderson-Higgs' mechanism by the spinon pairing.

Results
The paired hardcore boson model. We start from the paired hardcore boson model on the pyrochlore lattice, where the Hamiltonian is given as Here, b † i (b i ) creates (annihilates) one boson at the lattice site i, and n i ≡ b † i b i is the boson occupation number. This model differs from the usual hardcore boson model 9,26,[40][41][42] by having an extra boson pairing term. Previous theoretical works and numerical efforts on the hardcore boson model without the boson pairing have established the presence of the U (1) Coulomb liquid ground state that supports the gapless U (1) gauge photon and fractionalized excitations 40,43 . The main purpose of this work is to understand the role of this boson pairing on the phase diagram of the paired hardcore boson model.
This hardcore model has a strong physical motivation. This model is identical to the XYZ spin model via the standard mapping b i ≡ S − i , n i ≡ S z i + 1/2. The spin model was derived as a generic and realistic model that describes the interaction between the so-called "dipoleoctupole doublets" on the pyrochlore lattice 7,8,38 . The boson pairing naturally arises from the spin-orbit entanglement of the dipole-octupole doublets. In the end of this work, we further mention the relevance with the cold-atom systems that have been proposed 44,45 . Due to the boson pairing, the global U (1) symmetry is absent and the total boson particle number is not conserved, but the Hamiltonian remains invariant under a global Z 2 (or Ising) symmetry transformation with Throughout this work, we work on the regime with an average 1/2-boson filling. In the following, we first carry out the theoretical analysis and provide the physical understanding of the internal and emergent gauge structure and fractionalized excitations of this model, and then implement the large-scale QMC simulation to confirm the theoretical expectation.
The internal gauge structure and phase diagram.
Since the hardcore boson model without pairing is equivalent to the XXZ spin model and has been extensively Here 'expo' refers to 'exponentially', andñi ≡ ni − 1 2 . The ordered phase in the upper right region of the phase diagram in Figure 1 breaks the global Z2 (or Ising) symmetry. studied 9,40-42 , we briefly explain the ground state in the limit with t 2 = 0. When the hopping t 1 is greater than a critical value, the bosons are simply condensed and form a superfluid by breaking the global U (1) symmetry. In the opposite case when t 1 is less than a critical value, the system would form a U (1) Coulomb liquid with an emergent U (1) gauge structure and fractionalized excitations. Note the emergent U (1) gauge structure in the U (1) Coulomb liquid has nothing to do with the global U (1) symmetry of the model in the XXZ limit. Due to the emergent non-locality of the underlying U (1) gauge structure, the Coulomb liquid in the small t 1 regime is robust against any small and local perturbation such as the weak t 2 boson pairing.
The U (1) Coulomb liquid in the phase diagram can also be established from the limit with t 1 = 0. As we elaborate in the Supplementary materials, a sixth order degenerate perturbation theory in the t 2 pairing is needed to generate the three-boson hopping on the perimeter of the elementary hexagon of the pyrochlore lattice. It is this three-boson collective hopping that allows the system to fluctuate quantum mechanically within the extensively degenerate ground state manifold (or spin ice 9,46-48 manifold in the spin language) of the predominant boson interaction and lead to the U (1) Coulomb liquid. When both t 1 and t 2 are present and remain small, similar perturbative treatment again leads to U (1) Coulomb liquid. Therefore, we expect the U (1) Coulomb liquid to appear as the ground state when both t 1 and t 2 are reasonably smaller than V .
To establish the phase diagram, we first realize that the system favors a ferromagnetic order with S To reveal the connection between the Coulomb liquid and the ordered phases, we view the Coulomb liquid as the parent phase and implement the spinon-gauge construction 14,15 for the hardcore boson operators that is appropriate for the Coulomb liquid phase, where Φ † r (Φ r ) creates (annihilates) a spinon at the center (labeled by 'r') of the tetrahedron ('tet r '), and η r = ±1 for two sublattices of the diamond lattice formed by the tetrahedral centers. As we explain in details in Methods, the paired hardcore boson model becomes It is noticed that the boson pairing mediates the spinon interaction in the spinon-gauge formulation 14 . The spinon interaction may induce pairing between these fractionalized degrees of freedom and thus gap out the continuous part of the internal U (1) gauge field via an internal Anderson-Higgs' mechanism. Through the standard mean-field analysis, we do not actually find any pairing instability within the Coulomb liquid phase in the mean-field phase diagram (see Figure 1). This can be a mean-field artifact. Nevertheless, the mean-field theory does give a large region for U (1) Coulomb liquid in the phase diagram.
QMC algorithm. To examine the theoretical understanding, we perform the worm-type QMC algorithm 49,50 to simulate the model in Eq. (1). Since there are both boson hopping and pairing terms in the model, the typical worm QMC is no longer sufficient and a new update scheme is needed, which we outline here and more details can be found in Methods. We first express the partition function via Trotter product expansion in the imaginary time. We split the Hamiltonian into the free part c.] and the interaction part U = ij V n i n j − µ i n i , and expand the partition function with respect to the interaction part (using occupation basis that is denoted as |α ), where the grand canonical ensemble for the bosons is used by introducing the chemical potential to the bosons. The partition function is given as where k ≡ k 1 + k 2 . Under this representation, the configuration space of the partition function Z consists of all trajectories with closed worldlines (see Figure 2(a)), where "closed" refers to a periodic boundary condition with |α(0) = |α(β) . Due to the off-diagonal operator K, a boson can hop from one site to its neighbors via −t 1 b † i b j , or one pair of bosons can be created or annihilated at the same imaginary time and these two processes are dubbed hopping and pairing kinks, respectively. The numbers of such kinks are given by k 1 and k 2 in Eq. (4). To evaluate the dynamical properties, we further define a particular Green's function as As we show in Figure 2(b), G(i, τ I ; j, τ M ) introduces open trajectories that contain two worldline discontinuities "I" and "M", and (i, τ I ) and (j, τ M ) are the spatial, temporal locations of two worldline discontinuities. Shifting the discontinuities in space and time produces a series of trajectories. This is a crucial benefit of the worm-type algorithm that we can calculate the Green's function as efficiently as other thermodynamic quantities.
All closed and open trajectories constitute the total configuration space of the worm-type algorithm. Through three types of update procedures with certain probabilities, we can produce a Markov chain of different trajectories that walk in the total configuration space randomly. These procedures are classified as: (1) creation and annihilation of two worldline discontinuities,  in the following calculations. To determine the phase boundary between the disordered liquid phases and the ordered phase, we monitor the first order derivative of the free energy over the parameters t 1 and t 2 with We simulate these values by varying t 2 for fixed t 1 's with the system size N = 4 × 8 3 , β = (k B T ) −1 = 800, where we set V = 1 as the energy unit. The numerical phase diagram is presented in Figure 1(b). The transitions are strongly first order at small t 1 's and are consistent with the theoretical results in Figure 1(a). Moreover, as the system approaches the phase boundary near the horizontal axis, the transition becomes weakly first order like. In general, the phase boundary in Figure 1(b) is qualitatively consistent with the theoretical one.
To understand different phases, we probe the thermodynamic properties by measuring the specific heat and the entropy for the representative points in Figure 1(b). The results are depicted in Figure 4. For the U (1) Coulomb liquid in the pyrochlore ice context 9,10,48 , it is well-known that there exist double peaks in the heat capacity. The high temperature peak signals the entering into the spin ice manifold, while the low temperature peak arises from the quantum fluctuation that breaks the classical degeneracy of the spin ice manifold. Between the two peaks, there is an entropy plateau at the value of the Pauling entropy since the system is fluctuating within the ice manifold. Below the low temperature peak, the specific heat behaves as C v ∝ T 3 in the zero temperature limit due to the gapless U (1) gauge photon 15 . For the representative points 1,2 in Figure 1 the behavior of the specific heat is consistent with the U (1) Coulomb liquid (see Figure 4). This gapless excitation is the key signature of the emergent gauge dynamics, and is not related to any continuous symmetry breaking, especially since there is no symmetry breaking in the disordered regime and the (generic) model 7 does not even have a continuous symmetry. For the Z 2 liquid, all the excitations are fully gapped. Since the spinon pairing is expected to occur at very low energy scale, the double peaks in the heat capacity should persist except that we have an activated behavior of the heat capacity below the low temperature peak instead of the T 3 behavior for the U (1) Coulomb liquid. Inside the disordered regime of Figure 1(b), we find that the behaviors of "points 3,4,5" are consistent with a Z 2 liquid (see Figure 4 and Supplementary material). This result provides a thermodynamic evidence for the presence of a Z 2 liquid phase in the (orange) region of the disordered regime. More specifically, the thermodynamic gap, that is extracted from the heat capacity for the parameter point 4, is ∼ 0.018V . This is of the same order as the t 2 value, suggesting the possible physical origin of the Z 2 liquid state. As it was noted, the t 2 term renders an effective interaction between the (fractionalized) spinon quasiparticles. When one pair of spinons is condensed and individual spinon remains uncondensed, the U (1) Coulomb liquid would give way to the Z 2 liquid in a way similar to the superconducting pairing transition in a BCS superconductor. More physically, as t 1 /V increases inside the U (1) Coulomb liquid, the spinon gap monotonically decreases, and the interaction t 2 could lower the spinon pairing energy and overcome the reduced twospinon gap, leading to the Z 2 liquid state. In the Sup-plementary material, we provide more discussion about the detailed features of the specific heat in Figure 4 and discuss the possibility of charge density wave as an alternative explanation.
As listed in Table I, another important distinction between different quantum phases lies in the spatial dependence of correlation functions. Here we numerically measure the density-density and the boson-boson correlators that are defined as where r is the spatial separation between the lattice sites i and j. In the spin language, C n would correspond to the S z -S z correlator, while C b corresponds to the S + -S − correlator.
We first compare the correlations of the U (1) Coulomb liquid and those of the ordered state. As we depict in Figure 5, the boson density correlators for the parameter points 1,2 decay as a 1/r 4 power-law with the distance, and the boson-boson correlators decay exponentially. This is consistent with the prediction from the U (1) Coulomb liquid in which the density correlator at long distances and low energies 15,16 is mapped to the U (1) gauge photon modes 51 and the boson-boson correlator reflects the gapped fractionalized (spinon) quasiparticles 52,53 . In contrast, for the representative parameter point 4 inside the ordered state, the boson density correlator decays exponentially, and the boson-boson correlator saturates to a constant since the system develops the order in b by breaking the global Z 2 symmetry and simultaneously gives rise to a gap for the density correlator.
For the Z 2 liquid state, all the correlators should decay exponentially with the spatial separations. In our calculation, we find the boson-boson correlation does indeed decay exponentially. For the density correlators, despite the thermodynamic gap, we were unable to show more convincingly the exponetially decaying behavior due to the finite system size in our simulation and the tiny energy gap. To resolve this, one may need even larger system sizes to carry out the simulation in the future work.

Discussion
We discuss the physical realization of our spin or hardcore boson model. The solid-state realization has been proposed for the dipole-octupole doublets and studied in the previous works by one of us and collaborators 7,8,38 . Several Nd-based [29][30][31][32][33][34][35][36][37] and Sm-based 54 pyrochlore magnets 55 have been proposed to realize the dipole-octupole doublets, though most of them seem to support magnetic orders with mixed dipolar and octupolar components [29][30][31][32][33][34][35][36][37] . The known example of spin liquid candidate is the Ce-based pyrochlore Ce 2 Sn 2 O 7 where the Ce 3+ ion gives a dipole-octupole doublet 38,39 . Therefore, Ce 2 Sn 2 O 7 should be a good candidate to examine the spin liquid physics of the XYZ spin model. Beyond the solid state context, the cold atoms on optical lattices can be used to realize exotic models such as the paired hardcore boson model in this work. In a previous proposal, Ref. 45 has designed a ring exchange interaction for the bosonic gases via a Raman transition to "molecular" states on optical lattices to simulate the U (1) lattice gauge fields, where this Raman coupling has the form φ † b i b j and φ refers to the "molecular" state. More recently, the cold alkali atoms stored in optical lattices or magnetic trap arrays were proposed to realize a broad class of spin-1/2 models including the XYZ model by admixing van der Waals interaction between fine-structure split Rydberg states with laser light. Following these early proposals, we suggest two cold-atom setups to realize our paired hardcore boson models. In the first setup, we closely follow Ref. 45 and also propose a resonant coupling of the bosons via a Raman transition to a "molecular" two-particle state. Instead of choosing the original d-wave symmetry to simulate the ring exchange in Ref. 45, we propose a s-wave symmetry and condense the molecular states φ. Such a design naturally gives rise to an (uniform) hardcore boson pairing term φ † b i b j for a given lattice. For the second setup, one can directly make use of the known results and methods in Ref. 44 and extend to other lattices.
In summary, we have studied a paired hardcore boson model (or XYZ spin model in the spin language) on a pyrochlore lattice and found the broad existence of extoic quantum ground states. We make various suggestions for the experimental realizations in the solid-state and cold-atom contexts. Mean-field scheme. We describe the mean-field description in some details so that the underlying gauge structure and spinon-gauge interaction can be manifest. In the following we use the hardcore boson and the spin languages interchangeably. The degenerate classical spin ice configuration is equivalent to the the occupuation configuration of two bosons on each tetrahedron. We start with the physical meaning of the boson operators b i and b † i . From the perturbative analysis that is explained in the Supplementary materials for completeness, we learn that, the three-boson collective hopping becomes the "magnetic term" in the U (1) gauge theory Hamiltonian and formulation. Thus, b i and b † i would correspond to the vector U (1) gauge link of the U (1) quantum electrodynamics from this perspective. Physically, the perturbative calculation restricts us to the low-energy classical spin ice manifold and 'throws' away the high-energy excited states. Clearly, perturbative effective Hamiltonian doe not have information about the (spinon) matter field that carres U (1) gauge charges. Applying b † i breaks the ice rule on two neighboring tetrahedra that share the site i. Spinon excitations are created on the diamond lattice that is formed by the tetrahedral centers. Thus, b i and b † i carry two pieces of physical content, and the spinongauge contruction 14,15 clearly reflects this.

Methods
In the spinon-gauge formulation in the main text (that was originally introduced in Refs. 14 and 15), we have enlarged the physical Hilbert space. To return to the physical Hilbert space, a constraint was imposed in the main text 14,15 ). Since the Q r counts the spinon number, we further have [Φ r , Q r ] = Φ r and [Φ † r , Q r ] = −Φ † r . The spinon-gauge formulation of the microscopic Hamiltonian has been introduced in the main text, and is written here with more detailed position indices for more readability, where we set the gauge link A rr = 0 since we are dealing with t 1 > 0 and zero-flux sector for the spinons (see Supplementary materials), e µ (µ = 1, 2, 3, 4) refers to one of the four nearest-neighbor vectors on the diamond lattice, and there is a double counting of µ and ν. The mean-field results are obtained by systematically decoupling the spinon interaction into two-spinon terms with self-consistent mean-field conditions. These procedures are standard and follow closely with Ref. 14 that deals with a different model for non-Kramers doublets on the pyrochlore lattice.
Quantum Monte Carlo. We here give some details about our worm-type algorithm. Procedure (1) and (2) are same as the conventional one. Procedure (3) there is no such procedure in conventional Bose-Hubbard model. Figure 6 gives the schematic diagram of these procedures.
There is a note worth discussion here. In every wormtype algorithm there is an arbitrary value of ω G which defines the relative weight of closed and open space. The form of detailed balance equation for Procedure (1) is as follows: where N s is the number of segments and it is proportional to N β approximately. N is the number of lattice site and β is the reciprocal of temperature. We can see that if there is no ω G and N β is very large then the acceptance of changing a open trajectory to an closed one will be very small and the algorithm will be very inefficient. Here we choose ω G = N β according to the common choice which makes ω G /N s ∼ const.. More details can be found in Ref. 49.

Data availability
The data that support the findings of this study are available from the corresponding authors (G.C. and Y.J.D.) upon request.

References
Supplementary Materials for "Extended Coulomb liquid of paired hardcore boson model on a pyrochlore lattice" Perturbation theory. For the completeness, we provide a perturbative analysis and understanding of our paired hardcore boson model. In the well-known limit without boson pairing (i.e. t 2 = 0), the third-order degenerate perturbation within the spin ice manifold generates a three-boson collective hopping on the elementary hexagon of the pyochlore lattice that is given by 9 where t coll = 12t 3 1 /V 2 is positive for t 1 > 0, and 1, 2, 3, 4, 5, 6 are the six lattice sites on the perimeter of the hexagon. In the opposite case with t 1 = 0 and t 2 = 0, we need a six order perturbation theory. Within this low-energy manifold, one then expresses the hardcore bosons as b † i ∼ e iA rr where A rr is the U (1) vector gauge potential on the link connecting the centers of neighboring tetrahedra, and the effective Hamiltonian simply becomes 9 where a positive t coll favors a zero-flux sector with curlA = 0 for the spinons and * refers to the elementary hexagon on the diamond lattice formed by the tetrahedral centers.
FIG. 7. The fourth order perturbation. Here we apply two "t2"pairing processes and two t1 hoppings. The location of t1 or t2 indicates the bond that the t1 hopping or t2 pairing is applied.
When both t 1 = 0 and t 2 = 0, we can have several mixed contributions from the t 1 and t 2 processes. For instance, one could apply t 1 processes twice and t 2 processes twice could generate the three-boson collective hopping (see Figure 7). All these cases at all orders of perturbation series give positive contributions to the collective boson hopping of t coll and thus do not change the sign of t coll (or the ring exchange in the spin language). This justifies the choice of the zero-flux sector for the spinon hopping on the diamond lattice in the Method.
More supporting data. In this part, we list additional QMC results to support our conclusion that was made in the main text. In Figure 8, we plot the specific heats C v and entropy densities S/R of the points 3,5,6 in Figure 1.
At low temperatures C v decays exponentially for the points 3 and 5. For the point 6, it is a power-law decay. There are entropy plateaus at the value of Pauling entropy 1 2 ln( 3 2 ) in the plots of the entropy curves. This suggests that all these three parameter points experience the degenerate classical spin ice manifold during cooling. The energy densities of the parameter points 1-6 with decreasing temperature are represented in Figure 9. Numerically the energy densities of the parameter points 2,6 show power-law behaviors and it is exponentially decaying for the parameter points 3,4,5. As for the point 1 the simulation is more difficult so the data below ∼ O(10 −3 T /V ) were not so great and we were unable to fit it in the plot.
More discussion about the specific heat. We notice that the low-temperature peak of the specific heat in Figure 4(c) is quite sharp and much sharper than the ones in Figure 4(b) and/or Figure 8(c). We provide a ther- modynamic explanation for this phenomenon. We start from the entropy plateau at the value of the Pauling entropy at an intermediate temperature, below which the entropy would be gradually lost as we cool the system. For a gapped system that is expected for Figure 4(c), the entropy loss of the low temperature regime would be relatively small due to the energy gap. In contrast, the entropy losses of the low temperature regime for a gapless case in Figure 4(b) and Figure 8(c) would certainly be more. As a result, from the conservation of entropy, the remaining entropy loss would take place near the low temperature peak, and we must have a larger entropy loss with a higher peak at the low temperature in Figure 4(c) to compensate the large remaining entropy.
Unlike the Z 2 liquid in 2D, the Z 2 liquid in 3D supports a finite temperature thermal transition from the thermal proliferation of the line-like extended excitations. If our proposal of Z 2 liquid does apply to the narrow region between the Coulomb liquid and the ordered phase, we would expect a thermal phase transition. Although we cannot resolve this due to the system size and numerical difficulty, it is possible that the low-temperature specific heat peak in Figure 4(c) could be associated with the thermal phase transition.
The possibility of charge density wave. In the narrow region between the Coulomb liquid and the Z 2 symmetry breaking state, we found a gapped state. In the main text, we discuss the result from the perspective a gapped Z 2 liquid state, and indeed our results are consistent with the expectation for a Z 2 liquid. Moreover, as we have argued in the main text, the spinons of the would-be and/or nearby U (1) Coulomb liquid have a very small energy gap, and the strength of the boson pairing could simply overcome this gap and gain energies from the spinon pairing. Despite that this is a quite reasonable account of the numerical results, an alternative explanation may also apply to this gapped regime, and we simply describe here. Although we think it is not very likely due to energetic reason that we explain below, it is possible that, in the narrow region, the system develops a charge density wave (CDW) order for the hardcore bosons. For this to occur, we need to have further neighbor densitydensity interactions to overcome the three-boson collective hopping that is the dominant low-energy process in the gapped regime. This requires higher order perturbation than the third order and is suppressed. Nevertheless, the CDW, if present, supports gapped excitations and a finite temperature thermal transition. The way to distinguish a Z 2 liquid from a CDW is to measure the densitydensity correlator for a large system to direct detect the translation symmetry breaking.