Polar Molecules with Three-Body Interactions on the Honeycomb Lattice

We study the phase diagram of ultra-cold bosonic polar molecules loaded on a two-dimensional optical lattice of hexagonal symmetry controlled by external electric and microwave fields. Following a recent proposal in Nature Physics \textbf{3}, 726 (2007), such a system is described by an extended Bose-Hubbard model of hard-core bosons, that includes both extended two- and three-body repulsions. Using quantum Monte-Carlo simulations, exact finite cluster calculations and the tensor network renormalization group, we explore the rich phase diagram of this system, resulting from the strongly competing nature of the three-body repulsions on the honeycomb lattice. Already in the classical limit, they induce complex solid states with large unit cells and macroscopic ground state degeneracies at different fractional lattice fillings. For the quantum regime, we obtain effective descriptions of the various phases in terms of emerging valence bond crystal states and quantum dimer models. Furthermore, we access the experimentally relevant parameter regime, and determine the stability of the crystalline phases towards strong two-body interactions.


Introduction
The interplay of competing interactions and quantum fluctuations is known to allow for interesting phases to emerge in many-body quantum systems. This route towards novel phases of matter has been explored intensively in recent years, in particular in the field of low-dimensional quantum magnetism [1], and in the context of ultra-cold quantum gases on optical lattices, where one gains a high degree of control of the interaction strength [2]. The dominant inter-particle potentials in such systems are typically described by two-body interaction and exchange terms. It appears fruitful to explore also realistic set-ups of many-body quantum systems that are dominated -via engineered interaction potentials -by multi-body interaction terms of e.g. three-particle type. Indeed, recently, polar molecules have been proposed as promising candidates towards realizing such many-body systems. Driven by significant progress towards producing degenerate gases of polar molecules [3,4,5,6,7], various proposals have been put forward, how to drive ultra-cold polar molecules into regimes of strong manybody interaction effects [8,9,10,11,12]. In a recent work [11] an effective interaction potential was derived for polar molecules on an optical lattice in the presence of static electric and microwave fields. It was found, that upon appropriately tuning the external fields, the interactions between the polar molecules become characterized by extended strong two-as well as three-body interactions.
Here, we derive the ground state phase diagram for these systems with dominating three-body interaction using quantum Monte-Carlo simulations, exact finite cluster calculations and the tensor network renormalization group. We find the presence of solid structures for unconventional filling factors with unit cells much larger than the period of the underlying lattice, and a macroscopic ground state degeneracy in the classical limit. In the quantum regime, these degeneracies of the low energy sector are lifted giving rise to valence bond crystal states.
The interaction potential for polar molecules within an opitcal lattice in the parameter regime, where three-body and two-body interactions are present takes the form with V ij = V /r 6 ij and W ijk = W/(r 3 ij r 3 jk +perm.). Here, r ij denotes the spatial separation between particles on lattice sites i and j, and n i the local density at site i. The the two-body interaction V and the three-body interaction W can be tuned to be of similar strength, V W [11]. In the following, we consider in particular the case of bosonic polar molecules. In this case, the suppressed tunneling of a particle to an already occupied site needs to be accounted for by a hard-core constraint on the bosonic occupations [11]. Since the extended interactions decay rapidly with the inter-particle separations (c.f. Figure 1), we truncate the interactions after the leading terms, resulting thus in an extended Bose-Hubbard model with two-and three-body nearest neighbour interactions only, Here, b i and b † i denote boson annihilation and creation operators respectively, and the local density operator n i = b † i b i has eigenvalues 0 and 1 in the hard-core limit. Furthermore, t denotes the nearest-neighbor hopping matrix element, and µ is a chemical potential, allowing to control the filling (i.e. the density) of the system between n = 0 (empty) and n = 1 (full).
In previous works, models of hard-core bosons with extended two-and three-body interactions as effective models of ultra-cold polar molecules were studied for the case of a one-dimensional optical lattice [13] and the two-dimensional square lattice geometry [14]. In the one-dimensional case, an incompressible phase at a filling n = 2/3 was established for dominant three-body interactions, stabilizing both change-density wave (CDW) and bond-order wave (BOW) long-ranged correlations, apart from conventional CDW phases appearing at half-filling (n = 1/2) in the presence of two-body interactions. On the square lattice, several solid phases at fractional fillings as well as supersolid phases were found in a semi-classical approximation, some of which could also be verified by full numerical simulations.
Here, we extend such systematic explorations to the case of the honeycomb lattice, where nearest-neighbour three-body repulsions lead to characteristic effects of strong frustrations. In order to explore the physics of this system, we used a combination of quantum Monte Carlo (QMC) simulations based on the generalized directed loop algorithm within the stochastic series expansion (SSE) representation [15,16,17], as well as exact finite cluster calculations and the tensor network renormalization group approach [18] for the classical (t = 0) limit of the above Hamiltonian, as detailed below. In the following Section 2, we discuss the results of our calculations on the ground-state phase diagram in the absence of two-body interactions (i.e. for V = 0), thus focusing on the main aspects of three-body repulsions on the honeycomb lattice. We observe states with large degeneracies in the classical limit, and the emergence of complex valence bond crystal (VBC) phases, to be described in detail below. In Section 3, we discuss the behaviour of the system as we perturb it by a finite two-body repulsion V , and explore the experimentally relevant parameter regime where V W [11]. We examine the stability range of the VBC phases, and find that a cascade of incompressible solid phases steams from the competition between two-and three-body repulsion terms.

Three-Body Interactions
In this section, we focus on the case V = 0 in order to explore explicitly the effects of three-body repulsions on the honeycomb lattice. The ground state phase diagram in Figure 2 summarizes the results from our analysis. It exhibits at low values of t/W a variety of incompressible phases of different (unconventional) fillings n = 9/16, 5/8 (= 10/16), 2/3 and 3/4 (= 12/16). Due to the incompressible nature of these incompressible phases, they lead to finitely extended plateaus in the µ-dependence of the density n, as shown e.g. for fixed t/W = 0.3 in the inset of Figure 2. The actual nature of these phases and the quantum phase transitions between them, will be discussed below.
For larger values of t/W 0.4, the system is eventually driven by its kinetic energy from these solid phases via first-order quantum melting transitions into a uniform superfluid phase with a finite superfluid density ρ s . In the QMC simulations, the superfluid density is obtained as ρ s = w 2 /(βt) from measuring the bosonic winding number w fluctuations in the standard way [19] (here, β = 1/T (k B = 1) denotes the inverse temperature). The first-order nature of the melting transitions follows from pronounced jumps that are observed upon crossing the transition lines in both the density and the superfluid density. As a typical example, we show in Figure 3 the behavior of n and ρ s near the quantum melting transition between the n = 9/16 incompressible phase and the superfluid phase at µ/W = 1. The first-order nature of the quantum melting transitions is indeed expected, since (as shown below) the incompressible phases break the space group symmetry, whereas in the uniform superfluid U(1) symmetry breaking occurs at T = 0.

Quantum Monte Carlo
In order to perform the SSE QMC simulations for the current model, we employed a cluster decomposition of the Hamiltonian in terms of trimers of nearest-neighbor sites, such that each cluster carries one of the three-body interaction terms. The addition of the nearest neighbor two-body repulsions V (in Section 3), proceeds in this decomposition scheme as well, with each two-body term being shared by four such trimers. In the SSE directed loop construction, we thus used a doubly-linked list of 6-leg vertices. While the algorithm performed well for large values of the hopping t, we were not able to reach significantly below t/W ≈ 0.1 due to the dynamical freezing in the Monte Carlo configurations, once the competing diagonal interaction terms dominate  the Hamiltonian (which is a general issue of the SSE for Hamiltonians dominated by large diagonal terms). In terms of the simulation temperature T , we find that an inverse temperature β = 1/T = 20 − 50W (k B = 1) provided an optimal trade-off between the finite temperature incoherence and the algorithmic performance (the SSE algorithmic cost scales linearly in β). Furthermore, in order to be commensurate with all superstructures identified in the incompressible phases, the linear system size L is required to be an integer multiple of 6 (the total number of lattice sites being N = 2L 2 , as the honeycomb lattice contains two sites per unit cell, forming the two sublattices A and B). Within the above temperature regime, we were able to simulate systems up to L = 24, in some cases up to L = 36 in linear extend. In phases with large superstructures -described below -we find that the autocorrelation times of the bosonic structures increase such that the algorithm is not able to tunnel between different realizations of the ordering pattern even within several 10 6 QMC sweeps, but resides in one particular sector of the ground-state manifold (which however varies upon performing independent runs with different random-number streams for a given set of model parameters). During a first stage of the thermalization process, we annealed the system in a cyclic way by heating it up and cooling it down slowly such that it is able to relax globally into one of the equivalent ground-states. A fixed temperature thermalization was then performed in the second stage of the thermalization process. This annealing approach yield better performance than parallel tempering over an extended temperature ranges, since we are mainly interested in the ground-state phase diagram of the system, where β is large.
We furthermore employed quantum parallel tempering [20] (in µ or t) at fixed low temperatures in certain regions of the phase diagram, in particular in order to study the quantum phase transitions between neighboring VBC phases. We discuss these aspects in more detail in Section 2.3. Our special analysis of the model in the classical limit (i.e. for t = 0) is presented in Section 2.4. There, we also assess the applicability of the tensor network renormalization group approach for classical systems to the current three-body repulsive model. In the following Section 2.2, we discuss in detail the nature of the incompressible phases that appear in Figure 2, and how they are characterized from our numerical analysis.

Incompressible Phases
The n = 9/16 VBC Phase: The first density plateau that is encountered upon filling the lattice has a density n = 9/16 and extends between 0 ≤ µ/W ≤ 2 in the classical limit (t = 0). It has an potential energy per lattice site (equal to the internal energy for t = 0) of E (9/16) pot = −9/16µ, and corresponds to the closest packing of hard-core bosons on the honeycomb lattice without introducing any three-body repulsions. From geometrical considerations in the classical limit, this bosonic structure is obtained by covering the lattice with equilateral triangles with side length 4 √ 2a (where a is the distance between two lattice sites), each covering 16 lattice sites. These triangles are filled by nine bosons in a staggered (checkerboard) arrangement in order to obtain the overall filling 9/16. Neighboring triangles differ in the placement of the checkerboard pattern. This leads to domain walls along the edges of the triangles, where pairs of bosons reside. For V = 0, this does however not lead to a potential energy penalty. For an illustration of this structure, see the left panel of Figure 4. This structure allows for a denser packing of the particles (and thus a higher filling) than the overall checkerboard state of filling n = 1/2, without introducing any three-body energy terms. On every forth hexagon (such as the hexagons indicated in the left panel of Figure 4), six triangles share a common corner. In the classical limit, the energy of the system remains unchanged, if two particles change their positions along such a hexagon by an angle of 2π/6, as indicated in Figure 4. This local move results in a classical ground state degeneracy W = 3 N/32 . Thus, the ground state entropy is extensive, and the entropy per site is S/N = (ln W )/N = log(3)/32 ≈ 0.034 in the thermodynamic limit.
For finite values of t, the local moves allowed in the classical configurations lead to resonances on the hexagonal plaquettes, corresponding to second-order hopping processes of the bosons, effectively rotating a hexagon by an angle of 2π/3, as illustrated in the right panel of Figure 4. The system is able this way to gain kinetic energy from these tunneling processes. Such resonances also provide the dominant quantum fluctuations on this density plateau, and stabilize a VBC phase with a superstructure of checkerboarded triangles, linked by the resonating hexagons. Through an order-bydisorder effect, the quantum dynamics thus selects the ground-state to be a coherent superposition of the local resonance states.
We can indeed identify these peculiar features of the 9/16 plateau phase from the QMC simulations. In Figure 5, we show the local density n i , along with the kinetic energy density per bond K ij , where K ij = b † i b j + b † j b i and sites i and j belong to a nearest neighbor bond on the honeycomb lattice, for a representative point within the n = 9/16 phase. We find the local density is close to n i = 1 and 0, thus exhibiting only few fluctuations, in a staggered (checkerboard) pattern within triangular structures. On the other hand, the density around a subset of the hexagons -those, where the triangular checkerboard patterns meet -is within an intermediate range, and it is along the bonds of these hexagons, where the kinetic energy is mainly located. Both observations indicate a residual density dynamics in this incompressible phase, located along these hexagons. The resonant nature of the corresponding hopping events is reflected in the bond-bond correlation function K ij K kl , shown in Figure 6. We identify the main contribution from the hexagonal resonances, which lead to the enhanced bond-bond correlation between the reference bond (marked by an ellipse) and the bonds atop and below the reference bond in Figure 6. Furthermore, correlations of the same strength are visible between the reference bond and its next-nearest neighbour bonds to the left and right. These result from residual quantum fluctuations that lead to the finite kinetic energy distribution inside the checkerboarded triangular structures in Figure 5. We do not observe long ranged bond-bond correlations, thus the hexagonal resonances evolve essentially independently of each other.  The left panel of Figure 7 shows the density structure factor for the 9/16 VBC structure, S( q) = 1/N ij n( x i )n( x j ) e i q( x i − x j ) . Here, n a ( x i ) denotes the local density operator at position x i . A 6-fold structure is identified, with one of the equivalent peaks positioned at q 0 = (π, 0). This structure relates to the superstructure of equilateral triangles and thus provides a characteristic feature in the density distribution of this phase. The right panel of Figure 7 shows a finite size scaling analysis of QMC data for the structure factor S 9/16 = S( q 0 ) at this characteristic wavevector q 0 for a system within the 9/16 VBC phase. We find that S 9/16 /L 2 indeed extrapolates in the thermodynamic limit (N → ∞) to a finite value, verifying the presence of long-range density order in the bosonic structure of the 9/16 VBC phase.
The n = 5/8 VBC Phase: The next density plateau encountered upon further increasing the chemical potential appears between 2 < µ/W ≤ 5 in the classical limit, and has a filling of n = 5/8. The classical potential energy equals E (5/8) pot = −5/8µ+W/8 in this phase. As for the 9/16 plateau, we obtain a VBC structure in the quantum regime. However, its effective description is more involved.
In the classical limit, each valid configuration at this filling can be mapped to a lozenge tiling of the two-dimensional plane. Each lozenge is formed by eight lattice sites, and contains the boson pattern shown in Figure 8a, with a trimer-dimer pair such that one three-body vertex is introduced. A boson covering of filling n = 5/8 and energy E (5/8) pot results whenever the full area of the lattice is covered with such lozenges. However, one has to introduce a parity constraint in mapping to the boson configuration: coloring the lozenges as shown in Figure 8a, only sides of the lozenges of equal color are allowed to touch, in order not to introduce additional three-body repulsions, which would lead to a higher potential energy. Furthermore, it is well known that each lozenge tiling of the plane is dual to a hard-core dimer covering of a honeycomb lattice. This honeycomb lattice consists of hexagons that are a factor of two larger in linear extend than the hexagons of the underlying honeycomb lattice in the bosonic model. Different lozenge tilings along with the corresponding dimer coverings and the underlying boson patterns, are shown in Figure 8b and c. This mapping from bosonic configurations to dimer coverings allows to obtain the ground state entropy of the n = 5/8 phase in the classical limit from that of closed-packed hard-core dimer coverings on the honeycomb superlattice, which equals S = 0.108 [21] (the additional freedom of how the honeycomb superstructure is embedded onto the underlying lattice, and the two possible ways of coloring the lozenges lead to a further, non-extensive factor of 3 × 2 to the ground state degeneracy). Among the various lozenge tilings we find the dice lattice, shown in Figure 8b, which corresponds to a staggered dimer covering on the honeycomb superlattice. Other tilings contain a special group of three lozenges shown in Figure 8c. Rotating the bosons along the central hexagon in this structure, as shown in Figure 8c, leads to another allowed configuration, with the group of the three lozenges being rotated. In the dimer covering, this leads to a flip of the dimer pattern around the hexagon, as also seen in Figure 8c. In c), we show two configurations that contain a group of three lozenges, where the rotation of the bosons within the central hexagon leads to a flip between the two shown configurations. The notation for the two corresponding states in the quantum dimer model is shown in the bottom row.
For finite t, a third-order boson hopping process corresponds to this local flip in the dimer configuration on the hexagonal plaquettes of the honeycomb superlattice, as illustrated in Figure 8c. Using degenerate perturbation theory in t, this local dimer-flip dynamics is described by a quantum dimer model with solely kinetic terms, and favors the formation of plaquette resonances as shown in Figure 8c. Here, |α p and |β p denote the two states on the hexagonal plaquette p on the superlattice in Figure 8c involved in the plaquette flip process. For finite values of t, the system tries to maximize the number of plaquette resonances, and the ground-state degeneracy is partially lifted [21]. From the analysis of the quantum dimer model on the honeycomb lattice [21], we thus find that the ground state of the bosonic system at filling n = 5/8 corresponds to the plaquette VBC phase of the quantum dimer model with purely kinetic terms. This phase is characterized by resonances among hexagons, and we indeed observe a corresponding pattern in the QMC simulations: to exhibit this, we show in Figure 9 the local density and the kinetic energy density on the nearest neighbor bonds for a representative point within the n = 5/8 phase. The data is taken at t/W = 0.2, outside the asymptotic regime where the perturbative derivation of the effective quantum dimer model strictly holds. Nevertheless, we can identify in Figure 9 a pattern that shares the structure of the plaquette VBC phase of the quantum dimer model, as indicated by crosses on the most dominant hexagonal resonances.
The n = 2/3 VBC Phase: Further increasing the chemical potential, the classical model exhibits a first-order phase transition at µ/W = 5 from the n = 5/8 phase to a n = 3/4 density plateau, with potential energy E (3/4) pot = −3/4µ + 3/4W , which (as well as E (5/8) pot ) equals −3W at µ/W = 5. This point in the phase diagram is highly degenerate; in particular, we find among the degenerate ground states also states with a filling of n = 2/3 and a potential energy E (2/3) pot = −2/3µ + 1/3W (equal to −3W at µ/W = 5). One particular such state is shown in Figure 10. It consists of parallel strips of nearest neighbor boson pairs, separated by zig-zag chains of occupied sites. Other states of the same filling and energy can be obtained from this configuration upon allowing the bosons along the chains segments to move to a neighboring site. One such move is indicated in Figure 10. Such processes however are not independent of each other: For example, if the change indicated in Figure 10 was made, the two bosons located at the next-nearest neighbors along the chain of the originally occupied site are blocked to perform a similar move, since that would lead to additional three-body repulsion terms, and thus a higher potential energy. These changes in the configuration are thus blocking each other.
We find from the QMC simulation, that an extended phase of filling n = 2/3 gets selected via an order-by-disorder effect out of the degenerate ground state manifold in the classical limit at µ/W = 5: a new plateau of filling n = 2/3 emerges in the quantum phase diagram of Figure 2. It vanishes in the classical limit, while upon increasing t, it extends well into the regime µ/W < 5. In order to analyse the nature of this emerging phase, we show in Figure 11 representative QMC data of the local density and kinetic energy density within the n = 2/3 phase. One identifies a rigid backbone of parallel strips of boson pairs, where the local density takes on values close to unity. Furthermore, between these strips we find sites with a reduced local density, which are linked perpendicular to the stripes' direction by bonds with an enhanced kinetic energy density. Such behavior in both the density and the kinetic energy distribution is in accord with bond resonances, induced by the processes illustrated in Figure 10. The opening of the n = 2/3 plateau for finite t/W can thus be understood given the large local kinetic energy that the system gains, and which outweighs the penalty in potential energy. Since in the classical limit the kinetic energy contribution vanishes, the potential energy penalty leads to the disappearance of this phase. In Figure 12, we furthermore show the bond-bond correlations in the kinetic energy for the n = 2/3 phase. While the correlations decay quickly in the direction parallel to the strips due to the previously mentioned blocking effect, they sustain over a rather wide range perpendicular to the stripe direction. This underlines the robustness of this emerging n = 2/3 phase in the quantum regime. Figure 11. QMC data of the local density n i (shading) and the kinetic energy density along the nearest-neighbor bonds K ij (line thickness and shading) for bosons on the honeycomb lattice in the n = 2/3 VBC phase at V = 0, t/W = 0.3, µ/W = 5, and a system with L = 12 at β = 20.
The n = 3/4 Solid Phase: Increasing the chemical potential beyond µ/W = 5 , the classical system enters a phase of filling n = 3/4, that concludes into the fully occupied state for µ/W > 9 + 3t/W . A particular classical state of this filling and a potential energy of E  Figure 13. Additional states of the same filling and potential energy can be constructed by applying global moves along parallel lines throughout the system, as shown in Figure 13. Such global moves shift all bosons on two parallel lines along one sublattice. The classical ground state entropy of this phase can be calculated by counting the number of states W obtained by such moves. As the number of lines to perform such moves is proportional to the linear system size L, the entropy per sites S = (ln W )/(2L 2 ) ∝ 1/L scales to zero in the thermodynamic limit.
Because these degenerate ground states are connected not by local moves, but by global displacements, we do not find resonant structures for finite t. Due to the global nature of the moves that relate the various classical ground states, we expect the QMC algorithm to generate states that are members of this rigid manifold. Indeed, from the QMC simulations, we obtain density patterns that show characteristics of the deformed superlattice of fully occupied hexagons shown in Figure 13. A representative result from the QMC simulations is shown in Figure 14. The classical degeneracy thus appears not to be lifted in the quantum regime.

VBC-VBC Quantum Phase Transitions
After having described the VBC phases appearing in the presence of three-body repulsions on the honeycomb lattice, we next turn to the quantum phase transitions between these incompressible phases. An exploration of the transition regions is challenging for the QMC algorithm, because neighboring VBC states are rather close in energy, with competing potential and kinetic contributions of similar size. To illustrate this issue, we consider a numerical example of the relevant energy scales. The potential energy difference between the two phases at 9/16 and 5/8 filling is ∆E pot = −1/16µ + 1/8W , which equals −0.0625W near µ/W = 3. The kinetic energy difference between the two VBC phases for a given value of e.g. t/W = 0.3 is obtained from the QMC simulations as ∆E kin (t/W = 0.3) ≈ 0.025W , and is thus of similar size. We find that the competition between these two similar energy scales leads to an extended transition region between the VBC phases, that gives rise to an apparent continuous increase in the density, as seen e.g. in the inset of Figure 2. We explicitly checked that such behavior occurs also for temperatures as slow as T /W = 0.05 and T /W = 0.01, i.e. well below the above energy scales. Furthermore, we find strong algorithmic hysteresis effects upon varying the chemical potential at fixed t/W through the transition region within a single simulation. Figure 15 exemplifies this behavior: Starting within the n = 5/8 phase and decreasing µ, the density is stable until reaching µ/W ∼ 2.1, where it drops to n = 9/16 rapidly. Starting from the n = 9/16 phase and increasing µ, the filling of n = 5/8 is not established even well inside the n = 9/16 VBC region. This behavior exhibits a metastability of the n = 9/16 state; the algorithm cannot perform effectively the re-arrangements that are necessary in order to establish the structure of the n = 5/8 phase, starting from a configuration typical for the n = 9/16 phase.
We employed quantum parallel tempering in µ/W in order to assess, if with the help of replica-exchanges between neighboring values of µ/W , the algorithm overcomes such metastabilities. Results of these calculations are shown in Figure 15 for systems of different sizes for V = 0 and t/W = 0.3. We find that the obtained density n mimics the behavior of the hysteresis curve. In particular, we observe an increase of n within the range of µ/W , where the upper hysteresis curve drops. In addition, the compressibility κ shows a peak within this region, that increases with system size. In the QMC simulations, the compressibility κ = ∂n/∂µ = β( n 2 − n 2 ) is obtained in terms of the density fluctuations. The density curve flattens for the larger values of µ/W , and does not reach 5/8 even at µ/W = 3. On the other hand, n reaches the value 9/16 for µ/W > 2. Hence, the replicas are not able to establish the n = 5/8 VBC structure, even though they tunnel repeatedly throughout the extended parameter range. These observations are consistent with a first-order quantum phase transition between the neighboring VBC phases.
However, one might be concerned, that the numerical difficulties could hide an intermediate phase that separates the VBC phases. In one such scenario, an intermediate phase is not commensurate with our lattice sizes. When simulating lattices of varying sizes (L=12, 24, 25, 36), we did however not detect any new structures in the transition regime. Of course, we cannot exclude phases with very large unit cells from our finite size study. Another possibility would be a superfluid or even supersolid phase in the intermediate region, such as observed in a similar model on the square lattice [14]. However, we can rule out a finite superfluid density in the transition region: Using quantum parallel tempering in t/W or driving a superfluid system (from large t/W = 0.5) towards the transition region at low t, resulted in ρ s always becoming zero for values of t/W 0.35. A third possibility for an intermediate phase would be a VBC "emulsion" i.e., a metastable phase mixture with domain walls separating 9/16 and 5/8 VBC-like domains. While we observe in the QMC simulations bosonic structures that show features of both the 9/16 and the 5/8 phase, which would be expected from such an VBC "emulsion" scenario, these structures are also consistent with a first-order transition, given the algorithmic metastability. Thus, we consider a first-order VBC-VBC quantum phase transition the most conservative scenario consistent with the QMC data and the apparent difficulty of the QMC algorithm in this parameter regime. While we cannot determine the precise location of the quantum phase transition, we take the peak position of the compressibility as an estimate. This is reflected in the phase diagram in Figure 2: The points denote the region where the phase transition occurs according to our estimate. Since we are however not able to provide the exact position, but an interval, we denote this interval by an errorbar. We note, that the above discussion applies likewise to the other inter-plateau transitions.

Classical Limit (t = 0)
Finite Cluster Studies: In order to analyse the ground state phase diagram of the present model, it is useful to consider also the classical limit t = 0 of the quantum Hamiltonian. In the particular case V = 0, the model reduces to the most simple classical model on the honeycomb lattice with extended three-body interactions. We are not aware of any previous numerical or even exact results on this statistical physics model. While one can construct states that appear in the classical limit from analysing boson configurations from the QMC simulations, we wanted to check the implications for the classical limit using independent methods. For this purpose, we solved the classical model on a finite hexagonal cluster of linear system size L = 4 exactly. Doing so, we confirm the extends of the densities plateaus as described above, including the absence of a n = 2/3 plateau in the classical limit. However, this system size suffers from the particular problem, that among the various states of the 5/8 plateau, only the bosonic state corresponding to the dice lattice lozenge tiling (c.f. Figure 8b) matches onto this finite cluster. While this particular tiling provides a valid boson covering in the classical limit, this restriction shadows the huge degeneracy of this phase. Namely, on the L = 4 cluster, no configuration can be realized, that contains the group of three lozenges shown in Figure 8c. In order to be commensurate with such coverings, the linear system size must be an integer multiple of 6. All other phases however are commensurate with this system size, and we can exclude different boson patterns than those described above for up to 32 sites exactly.
In order to check on larger clusters, whether the classical predictions for the bosonic structures are correct, we tried to employ classical Monte Carlo simulations. However, they suffer from dynamical freezing at the low temperatures required to explore the plateau structure. In that respect, the QMC simulations perform more efficient, and we thus used the QMC simulations in order to generate configurations of the classical model (this is possible within the SSE framework since on each propagation level one obtains a configuration of the classical model). While we do of course not sample these classical configurations with the correct statistical weight, we can nevertheless check if the various classical configuration obtained this way are consistent with our effective descriptions. In particular, we verified that each classical configuration of density n = 5/8 and potential energy E pot = E (5/8) pot indeed complies with the construction given in Section 2.2 in terms of lozenge tilings. We did not observe any classical configuration with the right density and potential energy, that would violate this construction. We are thus confident, that our understanding of the classical phases is correct.
Tensor Network Renormalization Group: Recently, an interesting novel approach to study thermodynamic properties of classical statistical models has been proposed. It is based on a tensor network representation of the partition function, and evaluates it directly in the thermodynamic limit using a renormalization procedure in the tensor network decomposition. For details about this method, we refer to the original paper by M. Levin and C. P. Nave [18], as well as recent works applying this approach to the Isingmodel on the triangular [22] and the Shastry-Sutherland lattice [23]. Since this method performed well in these cases, we tried to apply it also to our model of nearest-neighbor three-body repulsions on the honeycomb lattice. In fact, in the original publication, the method was described explicitly for tensor network models on the honeycomb lattice. We found, that one can indeed express the partition function of our model in terms of a tensor network. For this purpose, one specifies two cyclically symmetric tensors  Table 1. Non-zero tensor elements T A ijk for links i, j, k around a site of sublattice A, and T B ijk for links i, j, k around a site of sublattice B, respectively. All other finite tensor elements are obtained from the shown ones via cyclic permutation of the indices i, j, k.
function of a tensor network model on the honeycomb lattice is then given as i.e. Z is obtained from the product of all the tensors, contracting pairs of indices for each bond of the honeycomb lattice. In order for Z to match the partition function of the bosonic model (in the classical limit t = 0), the tensors T A ijk and T B ijk can be chosen as given in Table 1, with the dimension D = 4. Here, we considered the general case, where both W and V are finite. One verifies easily, that with these particular tensors, Z indeed recovers the partition function of the classical particle model. Once a representation of the classical model in terms of a tensor network has been obtained, we can proceed to perform the renormalization procedure to this tensor network. Doing so, involves an approximation, since within each renormalization step the dimension of the renormalized tensor network is truncated to a fixed maximum dimension D max > D. D max is thus a regularization parameter of the algorithm.From the approximatively evaluated partition function, one obtains the free energy F = − T N ln (Z) in the thermodynamic limit, from with the density n is calculated by taking numerical derivatives of F with respect to µ. Doing so for sufficiently low temperatures T , we obtain the low-T phase diagram from resolving the density plateau structure. For sufficiently large D max , the numerical data eventually converges to the final result. Within the tensor network renormalization procedure, one performs a singular value decomposition of a (D max ) 2 × (D max ) 2 matrix of contracted tensors, from which the renormalized tensors are constructed. After each step the tensors have be to be normalized such that all tensor elements are smaller or equal to unity to avoid overflows. Since this procedure is iterated many times, one needs to implement all matrix computations using high-precision floating point arithmetics (e.g. the initial non-zero tensor elements differ in about 18 orders of magnitude from for µ/W = 2 and β = 20/W at V = 0.). In Ref. [22], quadruple precision was found appropriate for an Ising model. Using a customized version of LAPACK [24] with 128 bit reals (floating point precision of about 10 −34 ), we find that for the current model this precision was still not sufficient in the relevant part of the phase diagram. While the method performs for the case of purely two-body interactions (i.e. for W = 0), the tensor iteration procedure does not converge for β W , once W dominates and the chemical potential reaches beyond µ/W ≈ 1. (For β < W , the method converged and the results could be verified by comparing to Monte Carlo simulations. However, for such high temperatures, the plateau structure is thermally smoothed out, thus providing no information about the ground state properties, which we are after.) As an example, we show in Fig. 16 the density curves that we obtained at β = 10. In agreement with our previous analysis, we still find that the system enters a n = 9/16 plateau. However, we cannot explore the full phase diagram using this method. We trace the failure of the tensor network renormalization group procedure to the fact, that the tensors for larger values of µ/W span many orders of magnitude due to the exponential dependence of the Boltzmann factor in the tensor elements. Thus a broad range of values is required in order to account for the physics of the model. However, due to the accumulation of round-off errors in the computation of the renormalized tensor network (which includes a large number of additions and multiplications), one loses the required precision in a numerical implementation of the algorithm after a small number of iteration steps. This problem appears to be a generic drawback of the tensor network renormalization group method, which we suffer from because of the broad range of magnitudes that our system implies in the tensor network representation. For a model with two-body interactions only (i.e. at W = 0), the range of values is significantly reduced, and in this case we could indeed recover the (well known) behavior of the model, which is equivalent to the ferromagnetic Ising model on the honeycomb lattice. As an example, we show in Fig. 17 the filling n as a function of µ/V at W = 0 and β/W = 10, with a pronounced plateau at n = 1/2 emerging.

Two-body Interactions
Thus far, we considered mainly the case of purely three-body repulsions, and explored the phases of the bosonic model in that parameter regime. However, from the derivation of the extended Hubbard model for ultra-cold polar molecules in Ref. [11] it is clear, that two-body interactions will be at least of the same strength as the three-body terms (the claim in Ref. [14] , that the interactions are solely of three-body type thus appears in contrast to the results in Ref. [11]). Hence, it is important to assess the influence of two-body terms on the physics of such models. Starting with the nearest-neighbour three-body term W , the next important interaction term to be considered is the nearestneighbor two-body repulsion V . The phase diagram for hard-core bosons with solely nearest-neighbor two-body interactions (W = 0, V = 0) features a half-filled checkerboard solid for small values of t/V < 0.5, surrounded by a superfluid phase, without any supersolid phases present in the quantum phase diagram [25]. In the current setup, we recover this checkerboard solid at large values of V /W . An important question is, whether the other phases discussed in Section 2 are stable towards the relevant parameter regime V W , and if new phases appear from the competition between two-and three-body interactions.

Cascaded Transition:
To make predictions about the impact of two-body interactions, we first consider the n = 9/16 phase in the classical limit. It consists of triangles with an edge length corresponding to the size of four hexagons (c.f. Figure 4). No threebody vertices are present in this state, but the structure has neighboring bosons along the edges of the triangles. For finite V , these boson pairs result in a potential energy penalty, which tends to destabilize the structure.
One possibility would be, that there is a direct transition from the n = 9/16 phase to the 1/2 solid upon increasing V /W . However, we find that one can construct intermediate states with densities 1/2 < n < 9/16, that are energetically preferred for certain ranges of V /W . We obtain such state, upon generalizing the classical n = 9/16 state, where checkerboarded (half-filled) triangles are separated by domain walls: we now consider the size of these triangular domains to vary. Namely, we consider a honeycomb lattice and cover it with equilateral triangles of edge length x (in units of the size of one hexagon). Each triangle thus covers x 2 lattice sites ( Figure 20 shows QMC data for the local density that corresponds to a configuration with x = 12). A staggered filling with bosons yields a lattice filling of n △ (x) = x(x + 1)/2 − 1. The boson pairs along the boundaries of the triangles cost an energy of P △ = V (3x − 4)/2, and we obtain the potential energy per lattice site as Minimizing the energy with respect to x gives x = 4 for V = 0 (we thus recover the states at filling 9/16) and x → ∞ as V /W → ∞ (we converge to the half-filled checkerboard state). The number of resonances in a system with N sites equals N/(2x 2 ), resulting in a macroscopic ground state degeneracy of entropy S/N = ln(3)/(2x 2 ) for these solid phases.
In order to derive the ground-state phase diagram based on this construction, we minimize the energy for fixed µ/W and V /W as a function of x. For µ/W ≥ 2 other competing states are the n = 5/8 and for µ/W ≥ 5 the n = 3/4 phase described in Section 2. We thus optimise the energy among these various states in order to determine the phase boundaries. The phase diagram obtained this way is shown in the left panel of in Figure 18. We find that in addition to the previously established phases, the system exhibits a whole cascade of new solid structures (with x running from 4 to infinity) that appear upon increasing the value of V /W . As an example, we show in the right panel of Figure 18 the filling as a function of V /W for a fixed value of µ/W = 1. Starting from the n = 9/16 plateau, a cascade of transitions eventually leads for V /W > 0.35 to the staggered plateau of filling 1/2. This cascade of plateau phases forms an incomplete devil's staircase, induced by the competing nature of the three-and two-body repulsion terms. Numerical Results: For finite t > 0, we expect deviations from the classical cascade structure, since the number of hexagonal resonances that appear for finite t decreases as 1/x 2 with increasing x. The system could thus skip some of the higher-x plateaus because of the stabilization of the low-x phases by their larger kinetic energy. This can result in entering the staggered phase after only a finite number of intermediate plateaus. Besides, quantum fluctuations stabilize the n = 2/3 phase and we have to account for the relevance of this new phase on the overall phase diagram. However, it is difficult to resolve most of the new solid structures numerically because they are not commensurate with our lattice sizes. Moreover, the differences in energy and filling for states of neighboring values of x are small, and the plateaus rather narrow already in the classical limit. We nevertheless obtain evidence from the QMC simulations for (i) a non-direct transition to the staggered solid upon turning on a finite two-body repulsion V , and (ii) the stabilization of new VBC phases. For this purpose, Figure 19 shows QMC data for the V -dependence of the filling n for µ/W = 4, 5 and 7 respectively, obtained at t/W = 0.3. There one clearly identifies a plateau corresponding to the x = 12 structure (of filling n = 77/144), commensurate with our lattice sizes of L = 12 and 24. In Figure 20, we show a representative boson covering obtained by QMC for µ/W = 4 at V /W = 1.25, which is in perfect agreement with the x = 12 structure introduced above and the formation of hexagonal resonances as in the n = 5/8 phase. From Figure 19, we find that the n = 2/3 phase is stable for small finite values of V /W < 0.2, before a transition takes place towards the n = 5/8 density plateau. Similarly, we find from Figure 19, that the n = 3/4 phase remains stable for small values of V . At larger values  of V /W we clearly resolve the n = 2/3 and the n = 5/8 plateau. While we thus obtain evidence for a cascaded transition to the checkerboard solid, we were not in a position to fully resolve this transitions within our finite size QMC simulations. We hence did not attempt to comply a full phase diagram for finite V in the quantum case, but considered the above set of representative cuts through the parameter space. However, considering the relevant interaction range V W [11], we still find that in particular the n = 5/8 VBC phase can be realized for such realistic values of the ratio between three-and two-body interaction terms.

Conclusion
We studied a model of hard-core bosons with strong three-body repulsions on the honeycomb lattice using quantum Monte Carlo simulations, exact cluster analysis and the tensor-network renormalization group approach. The system's ground state phase diagram exhibits besides a superfluid region several complex valence bond crystal phases at fractional fillings 9/16, 5/8, 3/2 and 3/4. The obtained quantum phase transitions between neighboring valence bond crystal phases are consistent with firstorder transitions, given mild energy differences in both the potential and the kinetic energy sectors. With regard to a possible experimental realization based on cold polar molecules in appropriately tuned external electric and microwave fields, we included in addition to the three-body repulsions also nearest neighbour two-body interactions, which in a realistic set-up are at least of the same strength as the three-body interactions. Considering the competition between the different interaction terms, we obtained a cascade of intermediate incompressible plateaus as the two-body interaction strength increases. Furthermore, we find that the valence bond crystal phase of filling 5/8 with an effective low-energy description in terms of a quantum dimer model remains accessible well within the reachable parameter regime. A further step towards modeling cold polar molecules on the honeycomb lattice would be to take into account the full long-ranged nature of both the three-body and two-body interactions [11]. Treating the longer ranged contributions appropriately appears to require the use of different calculational techniques, since already the leading contributions to both interaction sectors provided a challenge for the quantum Monte Carlo approach. Based on our current results, one could expect that longer-ranged interactions stabilize additional solid phases with large unit cells, at least in the classical regime. For finite hopping strengths, residual entropies of the classical states would be lifted, eventually leading to the emergence of complex resonating structures similar to those described above. Whether other exotic states e.g. with topological order can indeed be stabilized in three-body extended Bose-Hubbard models [11], thus far appears open to future investigations.