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

In nature we are surrounded with examples of ordered phases at low temperatures—e.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–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–17] and experimental [18–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₁–J₂ model, there are still open debates on whether the non-magnetic state present in this model is a valence bond solid [24–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 (BM) or Bose liquid. The hallmark feature of a BM is the presence of a singularity in momentum space, known as a Bose surface [29–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 BM 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 BM 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 (CDW) structure factor, the ground-state fidelity metric and the condensate fraction. Next, in section 3, we discuss the identifying characteristics of the BM 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 BHH model: Bose–Einstein condensate (BEC)–CDW, BEC–BM and BM (other phase)–CDW. The main results are summarized in section 5.

2.1. Models

The model proposed in [34] to exhibit a BM phase is the spin-1/2 frustrated antiferromagnetic-XY model on the honeycomb lattice

$$\begin{equation} H = J_1 \sum_{\langle{ij}\rangle} ( S_i^+ S_j^- + {\rm H.c.} ) + J_2 \sum_{\langle\!\langle{ij}\rangle\!\rangle} ( S_i^+ S_j^- + {\rm H.c.}), \end{equation} \tag{ 1 }$$

where S^±_i is an operator that flips a spin on site i and J₁ (J₂) is the nearest-neighbor (next-nearest-neighbor) spin exchange. In this model, the next-nearest-neighbor coupling introduces frustration as long as J₂ > 0 (antiferromagnetism).

The Hamiltonian in (1) maps to a hard-core boson model (S⁺_i → b^†_i, S⁻_i → b _i, and J_i → t_i)

$$\begin{equation} H = t_1 \sum_{\langle{ij}\rangle} \left( b_i^\dagger b_j^{\phantom\dagger} + {\rm H.c.} \right) + t_2 \sum_{\langle\!\langle{ij}\rangle\!\rangle} \left( b_i^\dagger b_j^{\phantom\dagger} + {\rm H.c.} \right). \end{equation} \tag{ 2 }$$

Here b^†_i (b _i) is an operator that creates (annihilates) a hard-core boson on site i and t₁ (t₂) is the nearest-neighbor (next-nearest-neighbor) hopping amplitude. This Hamiltonian was shown to feature four phases: a simple BEC (a zero momentum (k = 0) condensate), a BM (a gapless spin liquid) and two fragmented BEC states. The BM was found to be the ground state of this model over the parameter range 0.210(8) ⩽ t₂/t₁ ⩽ 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 BHH Hamiltonian [36]

$$\begin{equation} H = t_1 \sum_{\langle{ij}\rangle} \left( b_i^\dagger b_j^{\phantom\dagger} + {\rm H.c.} \right) + t_2 \sum_{\langle\!\langle{ij}\rangle\!\rangle} \left( b_i^\dagger b_j^{\phantom\dagger} \mathrm{e}^{\mathrm{i} \phi_{ij}} +{\rm H.c.} \right) + V \sum_{\langle{ij}\rangle} n_i n_j, \end{equation}\noindent \tag{ 3 }$$

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-model (S⁺_i → b^†_i, S⁻_i → b _i, n_i → S^z_i + 1/2, t₁ → J₁, t₂ → J₂ and V →J_z)

$$\begin{equation} \fl H = J_1 \sum_{\langle{ij}\rangle} \left( S_i^+ S_j^- + {\rm H.c.} \right) + J_2 \sum_{\langle\!\langle{ij}\rangle\!\rangle} \left( S_i^+ S_j^- \mathrm{e}^{\mathrm{i} \phi_{ij}} + {\rm H.c.}\right) + J_z \sum_{\langle{ij}\rangle} \left( S_i^z + \frac{1}{2} \right) \left( S_j^z + \frac{1}{2} \right) . \end{equation}\noindent \tag{ 4 }$$

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 BM 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 sets our unit of energy, and we fix t₂ = 0.3 to focus on transitions from the phase identified in [34] as a BM phase. This model has two limiting cases: (i) for V →∞, the Ising regime, the ground state is a CDW and (ii) for V =0 and ϕ = π, the non-frustrated regime, the ground state is a simple zero-momentum BEC with non-zero superfluid density.

2.2. 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–46]. The fidelity metric is a dimensionless, intensive quantity and is defined as

$$\begin{equation} g = \frac{2}{N} \frac{1 - F(\lambda, \delta\lambda)}{(\delta \lambda)^2}, \end{equation} \tag{ 5 }$$

where N is the number of sites and the fidelity F(λ,δλ) is

$$\begin{equation} F(\lambda, \delta\lambda) = \langle{\Psi_0(\lambda) | \Psi_0(\lambda+\delta\lambda)}\rangle, \end{equation} \tag{ 6 }$$

where 〈Ψ₀(λ)〉 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 CDW insulator, where one of the two sublattices is occupied while the other one is empty. This state spontaneously breaks the sixfold rotational symmetry down to threefold 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

$$\begin{equation} S_{\rm CDW} = \frac{1}{N} \sum_{i,j} \langle{(n_i^a - n_i^b)(n_j^a - n_j^b)}\rangle, \end{equation} \tag{ 7 }$$

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 BEC, 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 ρ_ij = 〈b^†_ib _j〉,

$$\begin{equation} f_{\mathrm{c}} = \Lambda_1 / N_{\mathrm{b}}, \end{equation} \tag{ 8 }$$

where Λ₁ 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, Λ₁ ∼ 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 one-particle state [49, 50], and the BM phase. In the former case, some of the largest eigenvalues are O(N_b) and could even be degenerate. For the BM, 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

$$\begin{equation} n({\bf k}) = \langle{\alpha_{{\bf k}}^\dagger \alpha_{\bf k}^{\phantom\dagger}}\rangle + \langle{\beta_{{\bf k}}^\dagger \beta_{\bf k}^{\phantom\dagger}}\rangle, \end{equation} \tag{ 9 }$$

where $\alpha _{{\bf k}}^{\phantom \dagger } = \sum _{i\in A} \mathrm {e}^{\mathrm {i} {\bf k} \cdot {\bf r}_i} b_i^\dagger b_i^{\phantom \dagger }$ and $\beta _{\bf k}^{\phantom \dagger } = \sum _{i\in B} \mathrm {e}^{\mathrm {i} {\bf k} \cdot {\bf r}_i} b_i^\dagger b_i^{\phantom \dagger }$ 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]

$$\begin{equation} \langle{n({\bf k})}\rangle_{\theta_x,\theta_y} = \oint \mathrm{d}\theta_x \oint \mathrm{d}\theta_y \langle{n({\bf k},\theta_x,\theta_y)}\rangle, \end{equation} \tag{ 10 }$$

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 [34, 36].

Zoom In Zoom Out Reset image size

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 BM (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 BM is the absence of any order and a singularity in the momentum distribution n(k). In figures 2(a) and (b), we show n(k) for two values of t₂/t₁ that are typical for the BM phase. For this phase, n(k) features a t₂/t₁-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₂/t₁, 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.

Zoom In Zoom Out Reset image size

In this work, we present evidence that, in the (ϕ,V ) plane (see (3)), the BHH model for t₂/t₁ = 0.3 exhibits (at least) three phases at half-filling. For strong coupling V , the ground state is a CDW, while (at least) two possible ground states exist at weak-coupling. In the frustrated regime (ϕ ∼ 0) at V =0, the system favors a BM, 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

$$\begin{equation} L^{-\gamma / \nu} S_{\rm CDW} = f[(V-V_{\mathrm{c}}) L^{1/\nu}], \end{equation} \tag{ 11 }$$

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 three-dimensional (3D) Ising [53] and XY universality classes [54] 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).

Zoom In Zoom Out Reset image size

We can check the robustness of this result by considering the condensate fraction, which scales [55–57] according to

$$\begin{equation} L^y f_{\mathrm{c}} = g[(V-V_{\mathrm{c}}) L^{1/\nu}], \end{equation} \tag{ 12 }$$

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 first-order 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 BM 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 BM and the BEC state, but our results are consistent with a second order or a weakly first-order transition.

Zoom In Zoom Out Reset image size

It is now interesting to study how the Bose surface changes as ϕ increases and one transitions between the BM and the BEC phase. In figure 5, we show the momentum distribution function for four values of ϕ and fixed V =0. As seen in figure 5(b), the Bose surface reduces in size as ϕ departs from zero and continues to shrink until condensation occurs at k = Γ (see figures 5(c) and (d)) for ϕ > ϕ_c.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

The large value of g at V =0 may be an indicator of a transition away from the BM at V =0⁺. It is somewhat similar to the behavior of g in both the one-dimensional 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 BM 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 = Γ.

Zoom In Zoom Out Reset image size

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.0 and t₂ = 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 (figures 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 BM 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.

Zoom In Zoom Out Reset image size

In summary, we have studied the phase diagram of the hard-core BHH Hamiltonian, which has allowed us to probe the effect of perturbations on the BM 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 BM 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 BM, (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 non-trivial, and could be an example of deconfined criticality.

We have also found that the BM is destroyed upon increasing V , before the Heisenberg point for nearest-neighbor interactions V =2J₁ can be reached. The presence of a next-nearest-neighbor repulsion may change this and result in transitions to other exotic phases.

This research was supported by NSF through JQI-PFC (CNV and KS), ONR (CNV and MR) and US-ARO (VG). The authors thank L Balents, M P A Fisher, T C Lang, and Z Y Meng for useful discussions.