Mott Insulators in a Fully-Frustrated Bose Hubbard Model on the Honeycomb Lattice

We examine the effects of quantum fluctuations on a classical spin liquid state in the fully-frustrated honeycomb lattice Bose Hubbard model using quantum Monte Carlo simulations. Frustration is induced explicitly in the model by modulating the sign of the interaction spatially around each lattice hexagon. A superfluid to Mott insulating quantum phase transition can be induced by varying the relative strength of the classical interaction and quantum hopping. In the cases where the interaction has a regular spatial modulation, hopping promotes a phase transition to a symmetry-broken valence-bond solid state. When the interaction is forced to have no regular pattern, the Mott insulating phase is found to be featureless and gapped, making it an interesting candidate state for a quantum spin liquid arising in a Hamiltonian with only nearest-neighbor interactions.


Introduction
A concentrated effort is currently underway to build exotic states of quantum matter in ultracold atomic systems [1][2][3][4][5]. These atomic systems are destined to become "quantum simulators" of condensed matter Hamiltonians, allowing us to extend our reach into new models, phases, and phase transitions inaccessible to traditional condensed matter theory or experiment [4,6]. Although efforts to construct complex Bose or Fermi Hubbard Hamiltonians is still a fledgling effort, experimental progress is rapid, motivating theorists to continue to study increasingly realistic models that may someday be built and used in the search for exotic quantum phenomena.
Concomitant with this revolution, there has been a resurging interest in frustrated systems recently, due to the possibility that when frustration suppresses the ordering tendencies of a system, more subtle exotic phenomena can be observed. Of particular success would be the stabilization of a quantum spin liquid [7][8][9][10] -a featureless Mott insulating (MI) phase with emergent gauge symmetry that harbors fractional excitations -in a cold atom system. Spin liquids are projected to occur in several two-dimensional Bose Hubbard models [11,12], however the typical feature shared by these Hamiltonians is that the range of the kinetic hopping (tunneling) and interactions are multi-particle, long-range, or otherwise prohibitively complicated for realization in real experimental setups.
In this paper, we construct a featureless Mott insulating phase from a Bose-Hubbard Hamiltonian with only nearest-neighbor interactions. Motivated by the resonating valence-bond (RVB) [7,8] phase contained in the triangular lattice quantum dimer model [13], our model is equivalent to a spin-1/2 XXZ model constructed by adding quantum XY (hopping) perturbations on top of a fully-frustrated Ising model on the honeycomb lattice. We find that the hopping term lifts the degeneracy of the classical Ising model ground state, promoting various Valence Bond Solid (VBS) phases (with coexisting charge-density wave order), depending on the spatial distribution of the lattice frustration. Only when the spatial position of the frustration is made random, essentially breaking all lattice symmetries, does an incompressible (gapped) Mott insulator, a candidate quantum spin liquid phase, appear in the low-temperature phase diagram of the model.

Fully-Frustrated Honeycomb Lattice Bose Hubbard Model
In a recent breakthrough, Becker et. al. [14] have constructed a versatile experimental setup for loading ultracold atoms into a triangular optical lattice, demonstrating for the first time a superfluid-MI transition in this 2D system. Their experimental realization involves creating a periodic potential constructed from three laser beams in the XY plane. A rotation of the polarization of the three lasers in the XY plane configures the potential minima in the geometry of a honeycomb lattice. In the typical unfrustrated case, this bipartite lattice is expected to induce antiferromagnetic ordering for atoms [15]. However, Becker et. al.'s demonstration that the sign of the tunneling matrix element can be reversed may suggest the ability to induce frustrated or anisotropic tunneling of the atoms between lattice sites. It may therefore be possible in the future to introduce frustration into otherwise geometrically unfrustrated models. We explore this scenario by examining the possible phases of the simplest fully-frustrated honeycomb lattice Bose Hubbard model with nearest-neighbor interactions using large-scale QMC simulations.
We consider very generally the simplest model for hard-core bosons on the honeycomb lattice: employing only nearest-neighbor hopping (b † i and b i are the boson creation and annihilation operators), and density-density interactions (n i = 0 or 1). Here, the density-density interaction is written such as to allow for the exact mapping to the quantum spin-1/2 XXZ Hamiltonian, and n i → S z i + 1/2. Considering first the limit when t = 0, one has the classical Ising model, where frustration maybe be induced by ensuring that the sign of J ij is modulated in such a way that nearest-neighbor interactions can not be fully satisfied around a hexagonal plaquette. In the extreme limit where each hexagon on the lattice is frustrated (i.e. the model is "fully-frustrated") the low-temperature phase is an extensively degenerate manifold of equal-energy states -a type of classical spin liquid [16]. Each configuration in the degenerate manifold of states can be mapped to a close-packed hard-core dimer configuration on the dual triangular lattice [17]. It is well-known that the quantum version of the triangular lattice dimer model contains an RVB liquid in its phase diagram [13]. Although an exact mapping from the quantum dimer model to a quantum spin model is not available, one may consider the t = 0 classical fully-frustrated Ising model groundstate as a starting point, and ask whether quantum perturbations in the t = 0 model contains the spin counterpart to the RVB phase [7,8]; i.e. a Z2 quantum spin liquid phase [13].
Our Hamiltonian 1 has a distinct advantage over some other frustrated spin models since it may be simulated without the sign problem using Stochastic Series Expansion (SSE) quantum Monte Carlo [18][19][20][21] at finite temperature. Imposing the fully-frustrated constraint but keeping the magnitude of the interaction isotropic on the lattice, one can re-write the model where A ij takes the value of zero or one on each bond, subject to the constraint that the sum over each hexagon, A ij , is odd. When t = 0 the model is not fully gaugeinvariant, and modulating the sign of the hopping t would introduce a sign problem; we nonetheless refer to the choice of A ij subject to the constraint a gauge choice.

DTS
DRS Random Figure 1. The three gauges of study in this paper, the "Discrete Translationally Symmetric" (DTS) gauge, the "Discrete Rotationally Symmetric" (DRS) gauge and the "random" gauge. Sold bonds have A ij = 0, and dashed bonds have A ij = 1. Unit cells are marked in red, and there is no unit cell for the random gauge.
In the classical limit, the choice of the gauge A ij does not affect the nature of the groundstate manifold, due to the one-to-one mapping to the classical hard-core dimer tiling, which retains a residual entropy of S = 0.214 per site [17]. In this paper, we examine three different gauge choices (each with J positive and negative) to study the effect of the quantum hopping (t-term) on the classically degenerate manifold (see Figure 1). The first choice of gauge we refer to as the "Discrete Translationally Symmetric" (DTS) gauge where the pattern of bonds where A ij = 1 reduces the symmetry of the lattice such that only a discrete translational symmetry remains. We also study a gauge we call the "Discrete Rotationally Symmetric" (DRS) gauge, where the A ij variables leave the system with a nontrivial rotational symmetry and reduced translational symmetry. The third gauge is one in which we use a classical Monte Carlo to generate an unbiased random pattern of A ij that satisfies the fully-frustrated requirement A ij =1, and this is called the "random" gauge. Figure 1 illustrates each gauge with their unit cell outlined (where a unit cell exists).
We now examine the results of QMC simulations on the models defined in the previous section. We fix t = 1 in all results to follow. Note first that when J/t = 0 the sign of J on the bonds is irrelevant, and the model becomes an XY-superfluid for all gauge choices, characterized by a non-zero superfluid density (or spin stiffness) of ρ s = 0.292(1). This is the limit where SSE QMC with directed loops works most efficiently, hence we approach the transition to Mott insulating behavior by increasing J/t from zero in the simulations. We begin by examining the superfluid density (measured via the winding number [22]) at a fixed low temperature, T = 0.1t, as a function of increasing J/t. Results are illustrated in Figure 2 for the three gauges. We observe that quantum phase transitions between the superfluid and various insulating states are promoted by an interaction |J c | ≈ 2. The specific value of the critical J c is similar, but not strictly identical for the different gauges, or for the opposite signs of J in a fixed gauge.
Using the three different gauge choices, we now have access to six different insulating states (including ±J). In order to determine the ordering nature of the ground state, if any, we first examine q-space structure factors at a fixed J and temperature T . |J| = 3 was chosen to be well in the insulating phase, but of sufficiently small interaction to ensure algorithmic ergodicity. The T dependence of most observables was examined to find the approximate low-T convergence: T = 0.1t was chosen in most of the below data. Also, although the SSE QMC is in general very efficient at finding the ground state of a system, it is also capable of becoming stuck in a local minima in configuration space, similar to classical Monte Carlo algorithms. An example of such a local minima would be a screw dislocation if the ground state pattern were broken into a set of strips. The energy cost of such a defect would be finite and although translating the defect should be possible, fixing it requires passing through intermediate states of even higher energy. In order to overcome this problem we use the technique of thermal annealing. We start the simulation at a high temperature, run a fixed number of Monte Carlo steps (enough to equilibrate) and then lower the temperature until the desired temperature is reached. By using this method we generate low temperature configurations that tend to be free of defects. Figure 3 illustrates the density-density structure factor, calculated as where r ij is the vector connecting the unit cell containing site i to the unit cell containing site j on a honeycomb lattice and N is the number of unit cells, and the sum runs over all sites of a particular sublattice. We use the vectors between unit cells and a sublattice decomposition so that all the necessary information is contained in the familiar Brillouin zone of the triangular lattice. Figure 3 shows the density-density structure factors of both signs of J for the three gauges of interest for the case where n i and n j occur on the same sublattice. In the case of the DTS and DRS gauges, sharp peaks occur at certain wavevectors q i in the q-dependent structure factor. A peak will represent long-range order if S(q i )/N survives in the thermodynamic limit. In Figure 4, we examine the finite-size scaling of the various peaks, which demonstrates that in the case of the DTS and DRS gauges, there is indeed long-range order in the particle density. We also measure the structure factor of the bond-bond correlation function to search for VBS order. This structure factor is defined as where a labels the bond connecting sites i and j and r ab is the vector connecting the unit cells containing bonds a and b. A VBS phase will have Bragg peaks in this structure factor, which survive the thermodynamic limit as discussed above. Figure 5 shows the structure factor of the bonds for all cases studied in this paper. Crystalline solid-like structures in the sites and the bonds are visible when the system finds in a particular symmetry broken state at low temperatures; these are the respective VBS phases. Next, we generate real-space images of the lattice in order to help us visualize the various ground state particle and bond patterns. In order to visualize quantum order in the 2+1 dimensional QMC simulation cell, we essentially "project" the SSE basis configurations and operator list back into the 2D plane -effectively averaging the  Figure 4. Scaling of the peaks of the density-density structure factor for the four ordered cases and one disordered case ("random" gauge with J = 3). Bragg peaks are signified by the function approaching a constant in the limit of 1/L → 0, where L = N/2 is the length of the system. Fluctuations in the scaling of the DRS gauges reflects the fact that the system has multiple ground states with distinct ordering wave vectors that must be searched, while fluctuations in the random gauge are a result of the chance presence of small clusters that support short-range correlations. In the thermodynamic limit such clusters may result in local ordering, but do not give rise to long-range order.
simulation cell over all Monte Carlo time. One then can imagine plotting averaged basis variables by associating a color with the expectation value of each site: red for n i = 0, blue for n i = 1, and mixed (shades of purple) for 0 < n i < 1 (Figure 6). Similarly, we average the expectation value of the kinetic energy operator B a , defined in Equation 6, over the imaginary time expansion, and associate a thickness of the bond proportional to B a .
The corresponding images of each gauge are presented in Figure 6. Figures 6(a),6(b),6(c) and 6(d) show clear long-range order in the sites and the bonds, matching up with the structure factors presented in Figure 3 and Figure 5. The lack of order in the structure factors for the random gauges is now clarified by the real-space images. In the case where J = −3 the system appears to be made up of random domains of filled and empty sites, consistent with a lack of long range order. The case where J = 3 is the other extreme in which the system is uniformly featureless in the density, and disordered in the bonds, both again consistent with the absence of long range order.
In all four cases the average filling is half (corresponding to zero magnetization in the spin language). The low temperature boson compressibility (uniform spin susceptibility) gives information about the density fluctuations in the groundstate: if the susceptibility goes to zero, the phase is a gapped (Mott) insulating state. The temperature at which it turns off then gives us information about the gap to particle excitations. In Figure 7, susceptibility data is shown for the DTS and DRS gauges for both signs of J, and detailed data collected when J = 3 for the random gauge. For the J = 3 random gauge, it is clear that at low temperature, the susceptibility is zero. At high temperatures (T ≈ 0.1), the value of the susceptibility is higher than other gauge choices. This is likely in part because there is no long-range ordered pattern that needs to be broken to insert a particle. Indeed, the lack of ordering in the structure factors and real-space picture, combined with the fact that at low-temperature the phase is incompressible, suggests that the state is a featureless Mott insulating state -a candidate quantum spin liquid.

Quantum Order by Disorder
The results from QMC show a superfluid-MI phase transition in each of the cases studied, where the Mott insulator is a symmetry-broken crystalline state whenever the gauge choice is repeating, and a featureless state when the gauge is random. Keeping in mind that the t term of the Hamiltonian is off-diagonal in the density-basis choice, we can understand the crystalline phases in the limit of large J/t. Namely, treating t as a perturbation on top of the classical Ising model degenerate groundstate, one can ask which local t operations can lower the kinetic energy (maximize hopping) without costing energy proportional to J. In all cases we can see that the crystalline (VBS) states discovered by the QMC satisfy these two criteria. The simplest case comes from the DTS gauge where J = −3. If we look where the expectation of the hopping operator is largest, it occurs on bonds where energy is minimized when a particle is adjacent to a hole. The two hopping sites are also connected to two frustrated bonds before and after exchanging their positions, as shown in Figure 8. The real-space figure of the ground state shows a mixing of the two states in Figure 8, suggesting that such a superposition is also energetically favorable in the simulation.
In the case where J = −3 in the manifold of states that satisfy the classical part of the Hamiltonian, the hopping operator will only act on bonds that prefer holes next to particles. This is because swapping such sites is the only way to remain in the set of states that satisfies the classical part of the Hamiltonian. The only scenario in which the hopping operator can act on an attractive bond is if it is frustrated. Looking at b † i b j Figure 8. The left state is taken to the right state by the action of the b † i b j element of the Hamiltonian. If the left state satisfies the classical (J) part of the Hamiltonian, so will the right state. In this way, the superposition of these two (local) states reduces the energy of the quantum (t) Hamiltonian without increasing energy in the classical part of the Hamiltonian. Blue (red) circles represent particles (holes), while the green zigzags represent frustrated bonds. Figure 9. When the central bond is a frustrated attractive bond, all action of the hopping operator has a state that costs and energy of order J (denoted by an "X") as either the state before and/or after hopping. Blue (red) circles represent particles (holes), while the green zigzags represent frustrated bonds. Figure 8 again, all the bonds connected to the two swapped sites, except the bond between them, switch from frustrated to unfrustrated and vice versa when the sites are swapped. Figure 9 shows how there is no configuration possible that would allow hopping on an attractive bond to be a part of the ground state.
By eliminating the possibility that the attractive bonds fluctuate to lower the system's energy we find that only the repulsive bonds will, and the numerical simulations confirm this result. The ground state can now be described as the superposition of states that allows each repulsive bond to fluctuate as much as possible. If we extend this simple idea to all the gauges and signs of J, we can check to see if the highly fluctuating bonds in each case correspond to local configurations that satisfy configurations similar to Figure 8. Figure 10 shows the local structure around each of the highly fluctuating bonds from the DTS and DRS gauges, and in each case the local configuration satisfies what we expect from the perturbation theory approach.
To further cement our understanding of the fluctuations, we can use the information extracted from QMC measurement of the bond-bond correlation function, B a B b , to see if it matches these expectations. We'll use the simplest case to build our intuition, the DTS gauge with J = −3. In this case when looking at a layer of highly fluctuating bonds in their dominant configuration, any of them may hop to reduce energy. Once one of them hops however, neighbors on the same layer that were able to hop before are no longer able to. When we compare this with numerical results, we find this nearest neighbor behavior is reproduced. Figure 11 shows the real-space image of the ground state generated using the bond-bond correlation function. Notice how the thickness of the bonds next to the reference bond (red), corresponding to B a B b , is thinner for those adjacent compared to those two unit cells away. The graph in Figure 12 illustrates the thickness as a function of position along the chain, and clearly shows that adjacent bonds are anti-correlated when compared to well separated bonds.

Discussion
From our QMC data we have demonstrated that a superfluid-Mott insulator phase transition can exist in the nearest-neighbor Bose Hubbard model on the honeycomb lattice where the interactions are constrained to be fully-frustrating around each hexagon, subject to a gauge choice. We are able to classify the Mott insulating phases in each case where the gauge is a repeating pattern as a symmetry-broken crystal with coexisting charge-density wave and valence-bond order (i.e. a CDW-VBS [23]). Bragg peaks in the density-density and bond-bond correlations functions confirm the long- Figure 11. A real-space picture of the lattice generated using the bond-bond and density-density correlation function, with the red bond and black site as the reference bond. Notice the subtle oscillation in bond thickness for those bonds on the same layer as the reference. range order that we have illustrated using a real-space averaging of the QMC simulation cell. One could describe the Mott insulating phases as occurring by an order-by-disorder mechanism, where the degeneracy of the classical manifold of states for t = 0 is lifted by the quantum "perturbation" or hopping. We have intuitively described this orderby-disorder mechanism by a simple perturbation theory picture to first order in the hopping, subject to the constraint that each hopping process costs no energy in the (classical) interaction.
When the gauge choice of the interaction frustration is not a repeating pattern, we find that the order-by-disorder mechanism is destroyed. In this case, QMC finds a superfluid-Mott insulating transition, where the Mott insulator is a disordered phase, featureless on long length-scales. Both the q-space structure factors and the real-space simulation cell pictures confirm the groundstate configurations have no long range order. For J = −3 (predominantly attractive interactions), repeated simulations show that the configuration of the ground state consists of domains of filled or empty sites. In this way, the random gauge with J = −3 can be considered a disordered solid in the density and valence bond sectors. Since most of the hopping dynamics take place on the interface between domains, excitations of this gauge should be similar to those in the DTS gauge when J = −3.
The random gauge with J = 3 is the most interesting and unusual case, lacking order in the density and displaying a ground state where every site is exactly half filling on average. The density-density correlation function quickly decays to zero (within error bars) in a few lattice spacings. Like all the other gauges, the random gauge also has a gap to single particle excitations, shown by the susceptibility reducing to zero at low temperatures. Despite the gauge choice breaking all lattice symmetries, this phase is not described by a Bose-glass, since it is gapped with a vanishing local density Edwards-Anderson (EA) order parameter. One should also note that although the MI has no long-range order, it does have a somewhat broad distribution of bond strength (as evident in the differences in bond thickness in Figure 6(f)), which is induced locally due to the random pattern of frustration. This local disorder can be thought of as occurring "on top" of the RVB resonances. One may also imagine a different phase with no EA order parameter but disordered valence bonds, where spins are coupled in randomly frozen pairs (singlets) -a so-called Valence Bond Glass [24]. In contrast, spins in our phase are not coupled into singlets but free to resonate around the lattice, suggesting its description as a "glassy" RVB phase.
Thus, one can consider this MI is a good candidate for a quantum spin liquid groundstate -a state which is often defined as a T = 0 disordered phase that has an emergent gauge symmetry, associated with fractional excitations [9]. It is therefore interesting to consider briefly the excitations out of the groundstate. In the limit of t = 0, the only excitations in our model correspond to inserting a single particle and raising the energy of the system by J. In the limit of small t/J, different superpositions of classically allowed states will have different expectations of energy, but finding eigenstates other than the ground state is exceedingly difficult in SSE QMC. Since the total particle number commutes with the Hamiltonian, we can use it as a quantum number to label different states. If we can change the total filling of our ground state without violating the classical constraint, the new state we have can only be different in energy by some polynomial of t, or it may have a zero energy difference. The earlier discussion of hopping taken with attractive bonds leads us to one method of constructing such a state.
Let us imagine a situation similar to Figure 8, except with the central bond taken to be a satisfied attractive bond. In this scenario if we insert or remove two particles the number of frustrated bonds does not change, but we are taken to new eigenstate of the quantum Hamiltonian. Using QMC simulations with an applied chemical potential we were able to explore the excited state with two extra particles, compared to a half filled system. The evolution of the simulation suggests a large barrier between the half filled state and the state with two extra particles, as the state with one extra particle must have an energy penalty of order J, consistent with our earlier discussion. Repeated simulations from random initial configurations show us that states away from half filling must have some energy penalty, as the simulation always finds states with half filling at low temperatures as the lowest energy state. This implies that there is an energy difference between the ground state and that with two extra particles, and the difference between these states is not of order J, but some order of t. It is therefore interesting that the lowest-energy excitations in this groundstate do not appear to be the typical vison and spinon expected of a Z2 spin liquid [11]. Further characterization of the excitations in the model will be needed to make a possibles connection to a theoretical description of a quantum spin liquid phase.
In conclusion, we have studied a fully-frustrated Bose Hubbard model on the honeycomb lattice, where the frustration is induced by explicitly modulating the sign of the interaction around each hexagon to preclude a satisfied interaction. This choice of (classical) gauge can induce order-by-disorder when perturbed by the quantum hopping term, in the case when the gauge choice has a regular spatial pattern. In the case where the gauge choice is explicitly chosen to break all lattice symmetries, a disordered Mott insulating or RVB state is induced. Further work to look at excitations in this phase may be needed in order to determine its suitability as a quantum spin liquid state, i.e. a phase with fractionalized excitations and emergent gauge symmetry. Although more refinement of the model may be necessary, the identification of such candidate spin liquid states in Bose Hubbard Hamiltonians with only nearest-neighbor interactions is an important step towards the creation of such phases in cold atomic systems in optical lattices.