Quantum phases of hard-core bosons in a frustrated honeycomb lattice

Using exact diagonalization calculations, we investigate the ground-state phase diagram of the hard-core Bose-Hubbard-Haldane model on the honeycomb lattice. This allows us to probe the stability of the Bose-metal phase proposed in Varney et al (2011 Phys. Rev. Lett. 107 077201), against various changes in the originally studied Hamiltonian.


Introduction
In nature we are surrounded with examples of ordered phases at low temperaturese.g. crystalline solid structures, magnetically ordered materials, superfluid and superconducting states, etc. While it is straightforward to think of these ordered phases melting as the temperature is increased into the familiar classical liquid or gaseous states that are commonplace in every aspect of our lives, it has been a long-standing question as to whether 'quantum melting' at zero temperature can act similarly to thermal effects and prevent ordering. For a quantum spin or boson system, the resulting state of matter is known as a quantum spin liquid [1]. The interest in such a hypothetical spin liquid has remained strong for decades, most prominently due to the discovery of high temperature superconductivity [2,3].
Of critical importance is whether a two(or higher)-dimensional system can host a quantum spin liquid. At present, there exists a complete classification of quantum orders [4], which divides hypothetical spin liquids into several distinct classes. Some theoretical stability arguments have also been presented showing that there is no fundamental obstacle to the existence of quantum spin liquids [5]. Gapped spin liquid phases have been observed in dimer models [6][7][8], and also a family of special exactly-solvable toy models were discovered which can support gapped and gapless spin liquid phases [9]. Although these discoveries clearly demonstrated that a spin-liquid phase may appear in two (or higher) dimensions, at least in toy models, whether the same type of exotic phase can appear in a realistic spin system remains unclear.
Very recently, there has been much numerical [10][11][12][13][14][15][16][17] and experimental [18][19][20] evidence to suggest the existence of gapped spin liquids in models with SU(2) symmetry, but it is still unclear why these simple models can support such exotic phases. Of particular note are the numerical discoveries of a gapped spin liquid in the Heisenberg model on the kagomé lattice [11] and in the Hubbard model on a honeycomb lattice [10]. The existence of the latter is especially surprising and remains under debate [21]. The nature of this phase has been the subject of many works and it has been argued that next-nearest-neighbor exchange coupling is the mechanism responsible for the quantum spin liquid [13,14,22,23]. However, despite a number of numerical investigations into this J 1 -J 2 model, there are still open debates on whether the non-magnetic state present in this model is a valence bond solid [24][25][26][27][28] or a quantum spin liquid [12,15,16].
Gapless spin liquids, which may have low-lying fermionic spinon excitations that strongly resemble a Fermi-liquid state, have remained more elusive. Because of these excitations and because spin-1/2 models can be mapped onto hard-core boson models, some of these gapless spin liquids are often referred to as a Bose metal or Bose liquid. The hallmark feature of a Bose metal is the presence of a singularity in momentum space, known as a Bose surface [29][30][31][32][33][34]. However, unlike a Fermi liquid, where the Fermi wave vector depends solely on the density of the fermions, the Bose wave vector depends on the control parameters of the Hamiltonian and can vary continuously at fixed particle density.
In this paper, we follow up on the proposition of such a putative Bose metal phase in a simple hard-core boson (XY ) model on the honeycomb lattice [34], with an analysis of the stability of this phase against various changes in the Hamiltonian studied originally. First, we examine the dependence of the Bose wave vector on a (phase) parameter that makes the model transition between frustrated and non frustrated regimes. Next, we show that the phase identified as a Bose metal is stable to breaking of time-reversal symmetry and is present in the phase diagram of the hard-core Bose-Hubbard-Haldane (BHH) model, which features (at least) three phase transitions.
The remainder of this paper is structured as follows. In section 2, we define the model Hamiltonian, briefly discuss the Lanczos algorithm, and define the key observables used in this study: the charge-density wave structure factor, the ground state fidelity metric, and the condensate fraction. Next, in section 3, we discuss the identifying characteristics of the Bose metal phase in the context of the hard-core boson (XY ) model and show how the Bose wave vector evolves as the parameters are varied. In section 4, we discuss the three phase transitions that we can identify in the Bose-Hubbard-Haldane model: BEC-CDW, BEC-Bose metal, and Bose metal(other phase)-CDW. The main results are summarized in section 5.

Models
The model proposed in [34] to exhibit a Bose metal phase is the spin-1/2 frustrated antiferromagnetic-XY model on the honeycomb lattice where S ± i is an operator that flips a spin on site i and J 1 (J 2 ) is the nearestneighbor (next-nearest-neighbor) spin exchange. In this model, the next-nearestneighbor coupling introduces frustration as long as J 2 > 0 (antiferromagnetism).
The Hamiltonian in (1) maps to a hard-core boson model ( Here b † i (b i ) is an operator that creates (annihilates) a hard-core boson on site i and t 1 (t 2 ) is the nearest-neighbor (next-nearest-neighbor) hopping amplitude. This Hamiltonian was shown to feature four phases: a simple Bose-Einstein condensate (BEC) [a zero momentum (k = 0) condensate], a Bose metal (a gapless spin liquid), and two fragmented BEC states. The Bose metal (BM) was found to be the ground state of this model over the parameter range 0.210(8) ≤ t 2 /t 1 ≤ 0.356(9) [34].
To better understand the stability of the latter phase, we consider a strongly interacting variant of the Haldane model [35], the hard-core Bose-Hubbard-Haldane Hamiltonian [36] which reduces to (2) for φ ij = 0 and V = 0. Here, V describes a nearest-neighbor repulsion and the next-nearest neighbor hopping term has a complex phase φ ij = ±φ, which is positive for particles moving in the counter-clockwise direction around a honeycomb. Note that the Hamiltonian in (3) can be mapped to a modified XXZ- The complex phase φ plays two important roles. Firstly, for φ = nπ, time-reversal symmetry is explicitly broken. Therefore, we can use this control parameter to study the stability of the Bose metal phase against time-reversal symmetry breaking. Secondly, in the spin language, as we increase the value of φ from 0 to π, the sign for the next-nearest-neighbor spin-spin interaction is flipped from positive (φ = 0) to negative (φ = π), i.e., the next-nearest-neighbor spin exchange changes from antiferromagnetic to ferromagnetic. Since frustration in this model originates from the antiferromagnetic next-nearest-neighbor spin exchange, we can use φ to tune the system from a frustrated (φ = 0) to a non-frustrated (φ = π) regime, and thus it enables us to explore the role of frustration in stabilizing the BM phase.
In what follows, t 1 = 1 sets our unit of energy, and we fix t 2 = 0.3 to focus on transitions from the phase identified in [34] as a BM phase. This model has two limiting cases: (1) for V → ∞, the Ising regime, the ground state is a charge density wave (CDW) and (2) for V = 0 and φ = π, the non-frustrated regime, the ground state is a simple zero-momentum BEC with non-zero superfluid density (SF).

Method and Measurements
To determine the properties of the ground state of (3), we utilize a variant of the Lanczos method [37]. This technique provides a simple and unbiased way to determine the exact ground-state wave-function for interacting Hamiltonians. One limitation of the original algorithm is that the Lanczos vectors may lose orthogonality, resulting in spurious eigenvalues [38]. Orthogonality can be restored through reorthogonalization [39], which requires storing the Lanczos vectors. The storage needs can then be reduced utilizing a restarting algorithm, and the most successful techniques are the implicitly restarted [40,41] and the thick-restart Lanczos algorithms [42]. These two methods are equivalent for Hermitian eigenvalue problems, and here we utilize the thick-restart method for its simplicity in implementation.
A generic and unbiased way of determining the location of a quantum phase transition is related to the ground-state fidelity metric, g [36,[43][44][45][46]. The fidelity metric is a dimensionless, intensive quantity and is defined as where N is the number of sites and the fidelity F (λ, δλ) is where |Ψ 0 (λ) is the ground state of H(λ), and λ is the control parameter of the Hamiltonian. For strong repulsive interactions, the ground state of the BHH model is a chargedensity-wave (CDW) insulator, where one of the two sublattices is occupied while the other one is empty. This state spontaneously breaks the six-fold rotational symmetry down to three-fold but leaves the lattice translational symmetry intact. In addition, because of the diagonal character of the order established, the structure factor that describes this phase is maximal at zero momentum. Thus, we define the CDW structure factor S CDW as where n a i and n b i are the number operators on sublattice a and b in the ith unit cell, respectively.
Another possible ordered state is a Bose-Einstein condensate, where, in our model, bosons can condense into quantum states in which different momenta are populated. According to the Penrose-Onsager criterion [47], the condensate fraction can be computed by diagonalizing the one-particle density matrix where Λ 1 is the largest eigenvalue of ρ ij and N b is the total number of bosons. In a BEC, the condensate occupation scales with the total number of bosons as the system size is increased, which is equivalent to stating that ρ ij exhibits off-diagonal long-range order [48]. Consequently, in a simple BEC, Λ 1 ∼ O(N b ) while all other eigenvalues are O(1) [49]. Aside from a simple BEC, the eigenspectrum of the single-particle density matrix can signal fragmentation, where condensation occurs to more than one effective oneparticle state [49,50], and the Bose metal phase. In the former case, some of the largest eigenvalues are O(N b ) and could even be degenerate. For the Bose metal, however, all of the eigenvalues of ρ ij are ∼ O(1). Thus finite size scaling of f c can help pinpoint the presence or absence of condensation. Further understanding of the latter two phases can be gained by calculating the single-particle occupation at different momentum points where α k = i∈A e ik·r i b † i b i and β k = i∈B e ik·r i b † i b i are boson annihilation operators at momentum k for the A and B sublattices, respectively. In order to minimize finite-size effects and fully probe the Brillouin zone we average over 40 × 40 twisted boundary conditions [51,52], where θ α is the flux associated with the twisted boundary condition. For any finite-size calculation, there are a large number of clusters that one could study, each with slightly different symmetry properties. In this work, we focus solely on clusters that can be described by a parallelogram or "tilted rectangle". The clusters used in this study are illustrated in figure 1 and are discussed in more detail in [36] and [34].

XY Model
In a previous study [34], we reported that the phase diagram of the XY model (1) on a honeycomb lattice has three quantum phase transitions separating four distinct phases. The four phases are: (i) a BEC k = Γ (antiferromagnetism), (ii) a Bose metal (spin liquid), (iii) a BEC at k = M (a collinear spin wave), and (iv) a BEC at k = K (120 • order).
The key signature of a Bose metal is the absence of any order and a singularity in the momentum distribution n(k). In figure 2(a) and 2(b), we show n(k) for two values of t 2 /t 1 that are typical for the Bose metal phase. For this phase, n(k) features a t 2 /t 1 -dependent Bose surface, which, as a guide to the eye, we indicate by a dashed red line. In general, the Bose wave vector q B at which the maxima of n(k) occurs increases with increasing t 2 /t 1 , as shown in figure 2(c). We emphasize that the maxima in n(k) do not reflect Bose-Einstein condensation as they do not scale with system size.

BHH Model
In this work, we present evidence that, in the (φ, V ) plane [see (3)], the Bose-Hubbard-Haldane model for t 2 /t 1 = 0.3 exhibits (at least) three phases at half-filling. For strong coupling V , the ground state is a charge-density-wave (CDW), while (at least) two possible ground states exist at weak-coupling. In the frustrated regime (φ ∼ 0) at V = 0, the system favors a Bose-metal, while the unfrustrated regime (φ ∼ π) favors a BEC. Consequently, we find that there are three types of transitions: (i) BEC-CDW, (ii) BEC-BM, and (iii) BM(other phase)-CDW. Let us first consider the BEC-CDW transition driven by V at constant φ. In figure 3, we show the properties of the system for φ = π. In panel (a), we show the fidelity metric versus V , which has a smooth peak that grows with system size, indicative of a second-order phase transition (which would be unconventional in this case in which the system transitions between two ordered states) or a weakly first order transition. If the former is true, the structure factor would scale according to the rule: where N is the number of sites, L = N 1/2 is the linear dimension, and γ = ν(2 − η).
Because of our small lattice sizes, we cannot pinpoint the exact nature of this transition. For example, using a scaling analysis based on the 3D Ising [54] and XY universality classes [55] yield very similar results. In figure 3(b), we show the CDW structure factor scaled in accordance with the 3D XY universality class, resulting in V c = 3.71 (7). We can check the robustness of this result by considering the condensate fraction, which scales [56][57][58] according to where y = (d + z − 2 + η). This is illustrated in figure 3(c), resulting in V c = 3.73(3). This result is quite close to the one obtained using the structure factor. We stress, once again, that although this appears to be a second-order transition between two ordered states, finite-size limitations do not allow us to rule out the possibility of a weak firstorder transition or the existence of a small intermediate phase separating the BEC and CDW states. Next, we examine the properties of the model as one transitions from the Bose metal to the BEC state. In figure 4, we show the same quantities as in figure 3 (this time versus φ) for V = 0. The fidelity metric is plotted in figure 4(a), and peaks at approximately φ ∼ 0.88. In figure 4(b), we show the structure factor, which does not scale with finite size in either phase. Figure 4(c) depicts the scaled condensate fraction, yielding φ c = 0.84 ± 0.14, consistent with the peak in the fidelity metric. (Note that the 18A cluster experiences a level crossing for φ < φ c .) As in the previous case, we cannot make definite statements about the nature of the transition between the Bose metal and the BEC state, but our results are consistent with a second order or a weakly first order transition.
It is now interesting to study how the Bose surface changes as φ increases and one transitions between the Bose-metal and the BEC phase. In figure 5, we show the A third phase transition is expected as V is increased from zero and the BM phase is destroyed to give rise to the large V CDW phase. We illustrate this regime in figure 6 by plotting the fidelity metric (left panels) and CDW structure factor (right panels) for, (a) φ = 0, (b) π/12, and (c) π/6. In the left panels, for all three values of φ, one can see a sort of two-peak structure in the fidelity metric.  The large value of g at V = 0 may be an indicator of a transition away from the Bose metal at V = 0 + . It is somewhat similar to the behavior of g in both the onedimensional Hubbard model, where the Mott-metal-insulator phase transition occurs for the onsite repulsion U = 0 + [44], and the two-dimensional hole-doped t-J model [45], where d-wave superconductivity was seen to develop for a superconducting inducing perturbation with vanishing strength. Another possibility is that the Bose-Metal is stable for positive and small values of V , but a transition to another phase occurs when V < 0. The peak produced by such a transition would also explain the structure we see in g. We have also investigated this model with negative values of V , and found that large peaks are present in the fidelity metric for V < 0. The position of those peaks had a strong dependence on the cluster geometry. Hence, exactly what happens to the BM phase in the region V ∼ 0 is something that requires further studies, maybe with other techniques that allow access to larger system sizes and a better finite size scaling analysis.
For all clusters and values of φ depicted in the left panels in figure 6, one can also see a clear peak in the fidelity metric for finite values of V . This feature indicates the onset of CDW order. The structure factor S CDW , depicted in the right panels in figure 6, make apparent that for values of V beyond that peak, the CDW structure factor scales with system size.
In figure 7, we illustrate how the momentum distribution function changes in the presence of interactions at fixed φ = 0. n(k) is shown for V = 0 in panel (a) and as the interactions are increased in panels (b) and (c). In figure 7(b), one can see that the Bose-surface broadens as V increases and moves closer to k = Γ. Increasing the nearest-neighbor repulsion further, so that the system enters in the CDW phase [figure 7(c)], results in a momentum distribution function that peaked at k = Γ, albeit without condensation. Instead, the structure factor S(k) is sharply peaked at k = Γ.
A summary of our calculations for different values of V and φ is presented in figure 8 as the phase diagram of the hard-core BHH model at half-filling with t 1 = 1.0 and t 2 = 0.3. For φ > π/4, the boundary of the CDW phase was identified by the crossing points in the scaling of the structure factor (figure 3). The boundary of the BEC phase was identified by the crossing points in the scaling of the condensate fraction ( Figs. 3 and 4), and, for small values of V , also using the maximum of the fidelity metric for the largest systems sizes (figure 4). For φ < π/4, the CDW transition boundary was determined by the position of the maximum in the peak in the fidelity metric for the largest system sizes ( figure 6). Note that, in that regime, the Bose metal phase was found to be stable for V = 0. On the other hand, for V between 0 and the boundary of the CDW phase, the large value of the fidelity metric, as well as the behavior of several observables studied in that region, prevent us from making a clear statement about the nature of the ground state.

Conclusion
In summary, we have studied the phase diagram of the hard-core Bose-Hubbard-Haldane Hamiltonian, which has allowed us to probe the effect of perturbations on the Bose metal phase found in the frustrated XY model on a honeycomb lattice [34]. In particular, we explored the parameter dependence of the Bose wave vector and verified that the Bose metal is stable under the effects of time-reversal and chiral symmetry breaking. We identified three phases in the phase diagram of the BHH model; (I) a Bose metal, (II) a BEC, and (III) a CDW. The phase transitions between the different phases were identified utilizing the ground state fidelity metric, the CDW structure factor, the condensate fraction, and the momentum distribution. The BEC-CDW transition appears to be second order, although finite-size effects prevent us from ruling out the possibility of a weak first-order transition or the existence of an intermediate phase separating the BEC and the CDW states. If this transition is indeed a direct second order phase transition, the critical point would be highly nontrivial, and could be an example of deconfined criticality.
We have also found that the Bose metal is destroyed upon increasing V , before the Heisenberg point for nearest neighbor interactions V = 2J 1 can be reached. The presence of a next-nearest-neighbor repulsion may change this and result in transitions to other exotic phases.