# Coexistence of long-range orders in a Bose-Holstein model

###### Abstract

Exploring supersolidity in naturally occurring and artificially designed systems has been and will continue to be an area of immense interest. Here, we study how superfluid and charge-density-wave (CDW) states cooperate or compete in a minimal model for hard-core-bosons (HCBs) coupled locally to optical phonons: a two-dimensional Bose-Holstein model. Our study is restricted to the parameter regimes of strong HCB-phonon coupling and non-adiabaticity. We use Quantum Monte Carlo simulation (involving stochastic-series-expansion technique) to study phase transitions and to investigate whether we have homogeneous or phase-separated coexistence. The effective Hamiltonian involves, besides a nearest-neighbor hopping and a nearest-neighbor repulsion, sizeable double-hopping terms (obtained from second-order perturbation). At densities not far from half-filling, in the parameter regime where the double-hopping terms are non-negligible (negligible) compared to the nearest-neighbor hopping, we get checkerboard-supersolidity (phase separation) with CDW being characterized by ordering wavevector .

###### pacs:

71.38.-k, 67.80.kb, 71.45.Lr, 74.20.Mn## I Introduction

Whether diagonal long range orders [such as charge-density-wave (CDW) and spin-density-wave (SDW)] and off-diagonal long range orders [such as superconducting and superfluid (SF) states] can coexist homogeneously in correlated electronic systems is a central issue in condensed matter physics. Coexistence of superconductivity and CDW has been studied in three-dimensional systems blanton () (such as doped with or ), quasi-two-dimensional systems withers () (such as the dichalcogenide and ) and quasi-one-dimensional systems ong (); abbamonte () (such as the trichalcogenide and doped spin ladder ).

Supersolidity, defined as homogeneous coexistence of superfluidity and crystalline order, was theoretically proposed more than 4 decades ago leggett () and has been a subject of debate since then. The first experimental claim kim () of observing supersolidity in helium-4 further intensified the debate and enhanced the interest in understanding the phenomena. The general consensus chan (); prokofev1 () is that a supersolid (SS) is not realizable in a perfect hcp crystal of ; nevertheless, superflow can occur along vacancies that can collect near extended structural defects prokofev1 (); prokofev2 (). Thus, the occurrence of supersolidity in bulk solid helium-4 is being seriously questioned.

Cold-atom systems offer another opportunity for realization of supersolidity. Theoretically scalettar1 (); scalettar2 (); kedar1 (); kedar2 (); sheng (); sengupta (); sarma1 (), lattice bosons with various types of interactions in diverse geometries have yielded supersolidity. However, there has been no experimental creation of optical lattices with effective long-range interactions that produce supersolidity. Furthermore, experimental techniques to detect signatures of supersolidity also need to be developed for optical lattices sarma2 ().

There have been numerous studies of supersolidity involving hard-core-bosons (HCBs) scalettar1 (); scalettar2 (); kedar1 (); kedar2 (). A lattice model for quantum liquids, such as the interacting bosonic helium-4 at low temperatures, needs a hard-core constraint to account for the exclusion of occupation of more than one atom at each lattice point matsubara (); matsuda (). The behavior of the ground state and low temperature excitations of such systems are largely controlled by couplings with phonons. Furthermore, local Cooper pairs [comprising of two electrons (of opposite spin) at a site] can also be regarded as HCBs ng (). In Bismuthates, such HCBs couple to the breathing mode of the oxygen cage surrounding the Bismuth ions varma (); tvr ().

Here, in this article, we study a two-dimensional (2D) Bose-Holstein (BH) model for HCBs on a square lattice where they can hop to nearest-neighbor (NN) sites and experience the HCB-phonon interactions via a Holstein-type term. Previously, exact diagonalization calculations were done on this model sanjoy () for a small system (i.e., lattice) to study the resulting phase diagram. Here we use stochastic-series-expansion (SSE) based quantum Monte Carlo (QMC) technique to simulate large-size lattices so that various phases in the thermodynamic limit can be identified more clearly. Unlike the model, a SS is realized in our BH system due to non-negligible transport within the same sublattice. At densities not far from half-filling and at sufficiently large HCB-phonon couplings, phase coexistence occurs; furthermore, in the phase-coexistence region, the system tends to phase separate at stronger couplings.

Our paper is organized as follows: section II deals with a discussion of the BH Hamiltonian, its transformations and its mapping to the equivalent spin model. Section III covers a description of the numerical method and the observables employed to characterize the orderings. In section IV, we detail our results and the corresponding analysis, both with and without the presence of same-sublattice hopping terms. Lastly, in section V, we summarize our results and draw conclusions.

## Ii Formulation

The BH Hamiltonian is given by

(1) |

where and denote the annihilation operators for phonons and HCB particles, respectively, and is the number operator for HCBs at site . Furthermore, is the amplitude for hopping at NN sites denoted by and is the frequency of optical phonons mahan (). An effective Hamiltonian for the HCB particles is obtained by first transforming this BH Hamiltonian to the polaronic frame of reference (using the Lang-Firsov transformation) and then performing perturbation theory as detailed in Refs. sanjoy, ; sahinur, . Interestingly, second-order perturbation theory yields a two-step hopping which produces, in a 2D square lattice, both next-nearest-neighbor (NNN) and next-to-next-nearest neighbor (NNNN) hopping terms besides the usual NN hopping term in the Hamiltonian. Moreover, a NN repulsion also results from two-step virtual hopping back and forth.

So we get an effective Hamiltonian for HCB particles on a 2D square lattice sanjoy ():

(2) |

where and denote NNN and NNNN sites respectively; exp, , and ; and . In the regime , we can make the approximations and with the approximations becoming exact for . The small parameter is and is obtained from (see Ref. pankaj, for details); our perturbation analysis is done in the nonadiabatic regime ( ) and at strong coupling ().

Since all the hoppings are non-frustrated, a QMC simulation of the system does not suffer from the negative sign problem. We employ SSE technique in our QMC simulation and investigate the co-existence or competition of CDW and superfluidity in various regimes of the parameter space. Here we should mention that our model for HCBs is equivalent to an extended XXZ spin-1/2 Hamiltonian as shown below:

(3) |

where , , and stand for NN, NNN, and NNNN pairs respectively. Furthermore, the operators for the HCBs are related to those of spin-1/2 particles as: and . A comparison of the parameters in Eqs. (2) and (3) yields: , , , and . Now, the magnetization of the system can be tuned by using an external magnetic field; then, a term should be added to the Hamiltonian in Eq. (3) where the magnetic field is given in units of .

0.5 | 1.0 | 1.5 | 2.0 | 2.5 | 3.0 | |
---|---|---|---|---|---|---|

0.444 | 1.355 | 2.725 | 8.017 | 45.485 | 478.571 | |

0.415 | 0.970 | 0.960 | 0.647 | 0.395 | 0.255 |

Presence of hopping terms for HCBs indicates that superfluidity (i.e., spontaneous breaking of the global U(1) gauge symmetry) can exist in the system. On the other hand, a large interaction strength suggests the possibility of a CDW. Thus, our objective is to study the compatibility of these two long range orders. Now, these two orders can coexist either in a phase separated form or homogeneously as a SS. It should be pointed out that, a model on a square lattice does not show a thermodynamically stable SS phase for HCBs scalettar1 (). On the other hand, striped SS behavior is found away from half filling when NNN repulsion () is considered in the modelscalettar1 (); scal3 (); sengupta ().

## Iii Numerical calculations using SSE-QMC

We now give details of the SSE-QMC simulation of our model, or equivalently, our extended XXZ spin-1/2 model. Finding the phase diagram in the present problem requires exploring various limits of the parameters (including high anisotropy in our spin model). In our numerical computations, we used directed loop update for efficient sampling of the configurations sse1 (); sse2 (). The ground state properties are captured by simulating at low enough temperatures, i.e., with being the linear dimension of the square latticescal2 (). We employ since our calculations with yield the same values for the observables (within the error bars of the calculations). From calculations involving various large system sizes, we infer the results in the thermodynamic limit.

As can be seen from the expressions of the two-spin matrix elements for Heisenberg spin models used in various SSE-QMC studies sse1 (); sse2 (), a positive parameter is introduced to ensure the positivity of all the matrix elements (see Appendix A for details). This is necessary so that they can be treated as probabilities. The value of can also affect the autocorrelation time of the desired variable. We found that keeping the numerical value of equal to at least a quarter of the anisotropy parameter , particularly near the transition region, helps keep the data in various bins uncorrelated when each bin contains data from large number of Monte Carlo sweeps (i.e., at least 15,00,000); for the value of used in the XXZ model, see Ref. sse1, .

In this work we are concerned with the diagonal order parameter [i.e, the structure factor at the Neel ordering vector ] and the off-diagonal order parameter of the SF density . A general expression for (for our HCB system) is given as

(4) |

where denotes the ensemble average. Since (or in our spin model) are diagonal in the basis, a QMC average can be computed easily.

The SF density sandvik () is given by , where is the free energy in the presence of twisted boundary conditions with angle of twist . This is an off-diagonal order parameter. In a QMC calculation, the SF density along -direction is calculated using where and represent the total no. of Hamiltonian operators transporting spin in the positive and negative -directions, respectively sandvik ().

A few benchmarking comparisons of our SSE results with exact results for a model are shown in the supplementary material.

## Iv Analysis of results and technicalities

Our calculations for the 2D model aim at obtaining the ground state phase diagram of the system and also at extending (to the thermodynamic limit) the findings presented by a Lanczos study on a small cluster in Ref. sanjoy, . We consider in the present study.

To obtain the phase diagram for our system, the interplay between the diagonal and off-diagonal orders and the transition between them needs thorough investigation. As antiferromagnetic order breaks SU(2) symmetry whereas a SF phase breaks U(1) symmetry, a phase transition between them (according to Landau theory) cannot be of second-order type. Furthermore, a transition to a phase separated state occurs when the system undergoes a first-order transition.

We work in the grand canonical ensemble where there is no constraint on the HCB particle number (or the magnetization in the equivalent extended XXZ model). In Fig. 1, we show the variation of , , and with magnetic field (expressed in units of ) for g=1.75 for both our extended model and its XXZ version (i.e., ) in square lattices with and . The comparison shows that the results for different size lattices almost coincide. At small values of the magnetic field , the system manifests half-filling in the case of HCBs (or zero magnetization in the case of the equivalent spin model); the CDW phase is formed with maximum values of with SF density simultaneously assuming zero value. At large values of , before the system is at complete filling, decreases to zero while becomes finite manifesting the SF phase. In the intermediate magnetic field region, phase transition occurs with both the orders coexisting in our extended XXZ model. The plot of magnetization [or () where is the density] in Fig. 1(b), indicates a discrete jump during the transition in the case of the XXZ model. This is a typical signature for a first-order transition. Contrastingly, a continuous smooth increase in is observed for the extended XXZ model clearly ruling out the possibility of a phase separation. Calculations using a canonical ensemble shows an inhomogeneous phase coexistence due to phase separation (PS) sanjoy ().

In the region adjacent to in Fig. 1(a), where overlap of and are observed, the system displays a SS phase. Analysis of various large system sizes confirms the picture in Fig. 1(a) that there is indeed a homogeneous coexistence of CDW and SF long-range orders in the system.

Next, we will do a further detailed case-by-case study for the system with and without the effects of same-sublattice (i.e., NNN and NNNN) hoppings.

### iv.1 Considering only NN hopping ()

Here we study the case of in Eq. (2), , the bare XXZ model. Calculations are done on a lattice. We find that the system loses its CDW order at half-filling for values of the coupling below the critical value of corresponding to the Heisenberg point of the XXZ model. For smaller values of , superfluidity develops for all values of filling between 0 and 1. The phase diagram for the XXZ model is depicted in Fig. 2(a). Our calculated phase diagram of the XXZ model is compatible with Fig. 1, Fig. 2a and Fig. 2b of Ref. schmid, . At half filling (i.e., or ), the Heisenberg point denotes the boundary between CDW and SF phases.

Fig. 3(a) shows the jumps in magnetization as soon as we increase beyond . The magnitude of the jump increases as the coupling increases. This jump implies a first-order transition and indicates a phase separated coexistence of the CDW and SF phases. A histogram analysis can also capture the jump in values. In Fig. 3(b), plotted for , the histograms change when magnetic field values are varied. For instance, when the magnetic field is set at , even when large number of Monte Carlo sweeps were used in our QMC computation. On increasing the magnetic field to , the values start showing a double-peaked structure (with one peak at and another peak at ) which is indicative of phase separation [see also inset in Fig. 3(b)]. Then, at a slightly higher magnetic field of , all the values seem to be centered around a mean value close to 0.18. This manifests the first-order transition. All the discontinuous transitions are due to inhomogeneous coexistence of the CDW and SF phases.

### iv.2 Considering all hoppings

The phase diagram for our BH model [obtained from Eq. (2)] is depicted in Fig. 2(b). On including the effects of NNN and NNNN hoppings (i.e., ) in the model of Eq. (2), there can be a difference in the densities in the two sublattices owing to NN repulsion and same-sublattice hopping. As shown in Table. 1, at intermediate values of (i.e., ), NNN hopping is comparable to NN hopping; consequently, a SS state can occur. Whereas at larger values of (i.e., for ), NNN hopping is fairly smaller than NN hopping and we can expect the same behavior as in the (or the XXZ) model. On account of particle-hole symmetry in our model, the phase diagram is symmetric about half-filling.

We will now examine various features of the phase diagram of Fig. 2(b). The phase diagram for our model was obtained by identifying the transition regions and the nature of the various phases. The variation of magnetization with magnetic field in a lattice is shown in Fig. 6 in the transition region. The results for show that the magnetization increases gradually without any jump as the magnetic field is increased. Hence, in Fig. 6(a), for CDW and SF phases coexist homogeneously resulting in a SS state. Here we must mention that Fig. 1(a) (depicting simultaneous coexistence of CDW and SF phases through non-zero values of and ) corroborates this conclusion.

As a matter of fact, large anisotropy [i.e., large values of in Eq. (2)] in the model requires large simulation time. Thus, as we increase the value of and thereby increase the value of or (see Table 1), the numerical calculation suffers from appreciable slowing down and with our computational constraints we cannot study for large values of (i.e., ).

We can set a cut-off for the anisotropy parameter above which the essential physics for the system does not change much; for values of above the cut-off, we expect that NN occupancy is in effect projected out. More precisely, we find that we can deal with large values of (i.e., ) by setting and still get the correct behavior of the observables thereby saving computational time. Fig. 4 shows versus and against for ; results for different values of are compared in order to decide on a cut-off value . In the large anisotropic limit, a change in the anisotropy parameter can be shown to produce an additional effective field of (at the mean-field level). So in Fig. 4(a), we shifted the scale accordingly in an attempt to make the plots coincide. The good agreement between the two cases for and gives us the freedom to use a cut-off of at large values of (i.e., ).

With the above simplification, we do the SSE-QMC simulation for larger values of (such as ) and observe phase-separated phases; we explain this based on Fig. 5 plotted at . At values of magnetic field , a single peaked structure occurs. On increasing h, at , a double-peaked structure results showing simultaneous existence of two phases (with magnetizations centered at and ). A further small increase to , leads to again a single peak (centered around ) signalling that a discontinuous phase transition has occurred. Thus a phase separated state at is clearly captured in Fig. 5(a). In Fig. 5(b), we show the evolution of the SF density and as is varied and capture the first-order transition at .

250289 | 96613 | |

658170 | 174995 | |

208295 | 9835 | |

960 | 113 |

238406 | 82608 | |

856138 | 666353 | |

161336 | 8646 | |

4847 | 4778 |

The results for the magnetization versus magnetic field in the transition region for values of and are shown in Figs. 6(b), (c), and (d). There is a gradual increase in the quantum of the magnetization jump as the coupling is increased. We should also mention here that, particularly in the calculations involving large , we wanted to ensure that the numerical data in adjacent bins are not correlated. To this end, we compute the autocorrelation time defined as sse1 ()

where

(5) |

with and being Monte Carlo times defined in units of Monte Carlo sweeps (MCS).

Typically, we see that for moderately high values of , can restrict from attaining very large values. Thus, for , we use a large bin size (i.e., 15,00,000 MCS) in our simulations and keep well within the bin size in order to produce meaningful results from our simulation. However, for larger , even choosing cannot keep autocorrelation times sufficiently smaller than such large bin sizes. So, in those cases, we use larger values (i.e., =6 and 8 for =15 and 20, respectively) and even larger bin sizes (i.e., 22,00,000 MCS). The values of the autocorrelation times, for and at various fields close to the transition, are shown in Table 2. It should also be noted that we cannot take too large, as a calculation with larger than does not produce meaningful results.

Lastly, we mention that our phase diagrams, both for the XXZ model and our extension of it, are similar (though not identical) to what were obtained using Lanczos method in Ref. sanjoy, .

## V Summary

In this work, we studied the effective Hamiltonian of a Bose-Holstein model using the SSE-QMC technique and obtained the ground state phase diagram. We found that supersolidity is realized at intermediate couplings; whereas, at large couplings the system phase separates because the double-hopping terms (that produce transport in the same sublattice) are not dominant.

Our results on a large lattice represent well the system behavior in the thermodynamic limit, as demonstrated in Fig. 1 using different finite-size calculations. Our results are only qualitatively similar to those obtained earlier using modified Lanczos technique on a much smaller cluster sanjoy ().

We overcame computational difficulties for large repulsive interactions by devising a cutoff repulsive strength; above the cutoff, the system properties (as shown in Fig. 4) become essentially independent of the strength of repulsion (because repulsive interactions project out nearest-neighbor occupation). To mimic the results of statistically independent configurations, we considered sufficiently large values and kept the auto-correlations within acceptable limits.

Our work is an exercise in SSE-QMC study of a simple but important model which has significance in various fields. We hope that the present results will stimulate further investigations in allied areas such as frustrated quantum magnets in various geometries and at various magnetizations; coherence dynamics of excitons/spins in the presence of phonon environments (pertinent to quantum computation and artificial light harvesting); dimer formation and dimer correlations in Hubbard-Holstein model murakami (), HCBs coupled to multimode phonons zoller (), etc.

## Appendix A SSE bond operators

In SSE-QMC study of Heisenberg spin systems, the Hamiltonian is written as a bond Hamiltonian. Particularly, in our case we write where denote the NN, NNN, and NNNN bonds in our spin model, respectively. Each of such consists of the diagonal () and the off-diagonal ) parts and is given as with expressions

(6) |

where , , , and with the coordination number . In our model, a two-spin matrix element of any of these operators can never become negative.

## References

- (1) S. H. Blanton, R. T. Collins, K. H. Kelleher, L. D. Rotter, Z. Schlesinger, D. G. Hinks, and Y. Zheng, Infrared study of from charge-density-wave insulator to superconductor, Phys. Rev. B 47, 996 (1993); A. M. Gabovich, A. I. Voitenko, and M. Ausloos, Charge- and spin-density waves in existing superconductors: competition between Cooper pairing and Peierls or excitonic instabilities, Phys. Rep. 367, 583 (2002).
- (2) For a review, see R. L. Withers and J. A. Wilson, An examination of the formation and characteristics of charge-density waves in inorganic materials with special reference to the two- and one-dimensional transition-metal chalcogenides, J. Phys. C 19, 4809 (1986).
- (3) W. W. Fuller, P. M. Chaikin, and N. P. Ong, Superconductivity and charge-density waves in Ta- and Ti-doped , Phys. Rev. B 24, 1333 (1981).
- (4) A. Rusydi, W. Ku, B. Schulz, R. Rauer, I. Mahns, D. Qi, X. Gao, A. T. S. Wee, P. Abbamonte, H. Eisaki, Y. Fujimaki, S. Uchida, and M. Rübhausen, Experimental Observation of the Crystallization of a Paired Holon State, Phys. Rev. Lett. 105, 026402 (2010); P. Abbamonte, G. Blumberg, A. Rusydi, A. Gozar, P. G. Evans, T. Siegrist, L. Venema, H. Eisaki, E. D. Isaacs, and G. A. Sawatzky, Crystallization of charge holes in the spin ladder of , Nature (London) 431, 1078 (2004).
- (5) A. J. Leggett, Can a Solid Be ”Superfluid”?, Phys. Rev. Lett. 25, 1543 (1970); A. Andreev and I. Lifshits, Quantum theory of defects in crystals, Zh. Eksp. Teor. Fiz. 56, 2057 (1969); G. Chester, Speculations on Bose-Einstein Condensation and Quantum Crystals, Phys. Rev. A 2, 256 (1970).
- (6) E. Kim and M. H. W. Chan, Probable observation of a supersolid helium phase, Nature (London) 427, 225 (2004); Observation of Superflow in Solid Helium, Science 305, 1941 (2004).
- (7) D. Y. Kim and M. H. W. Chan, Absence of Supersolidity in Solid Helium in Porous Vycor Glass, Phys. Rev. Lett. 109, 155301 (2012).
- (8) M. Boninsegni and N. V. Prokofev, Colloquium: Supersolids: What and where are they?, Rev. Mod. Phys. 84, 759 (2012).
- (9) M. Boninsegni, A. B. Kuklov, L. Pollet, N. V. Prokofâev, B. V. Svistunov, and M. Troyer, Luttinger Liquid in the Core of a Screw Dislocation in Helium-4, Phys. Rev. Lett. 99, 035301 (2007).
- (10) G. G. Batrouni and R. T. Scalettar, Phase Separation in Supersolids, Phys. Rev. Lett. 84, 1599 (2000).
- (11) F. Hebert, G. G. Batrouni, R. T. Scalettar, G. Schmid, M. Troyer, and A. Dorneich, Quantum phase transitions in the two-dimensional hardcore boson model, Phys. Rev. B 65, 014513 (2001).
- (12) D. Heidarian and K. Damle, Persistent Supersolid Phase of Hard-Core Bosons on the Triangular Lattice, Phys. Rev. Lett. 95, 127206 (2005).
- (13) A. Sen, P. Dutt, K. Damle, and R. Moessner, Variational Wave-Function Study of the Triangular Lattice Supersolid, Phys. Rev. Lett. 100, 147204 (2008).
- (14) R. G. Melko, A. Paramekanti, A. A. Burkov, A. Vishwanath, D. N. Sheng, and L. Balents, Supersolid Order from Disorder: Hard-Core Bosons on the Triangular Lattice, Phys. Rev. Lett. 95, 127207 (2005).
- (15) P. Sengupta, L. P. Pryadko, F. Alet, M. Troyer, and G. Schmid, Supersolids versus Phase Separation in Two-Dimensional Lattice Bosons, Phys. Rev. Lett. 94, 207202 (2005).
- (16) V. W. Scarola and S. Das Sarma, Quantum Phases of the Extended Bose-Hubbard Hamiltonian: Possibility of a Supersolid State of Cold Atoms in Optical Lattices, Phys. Rev. Lett. 95, 033003 (2005).
- (17) V. W. Scarola, E. Demler, and S. Das Sarma, Searching for a supersolid in cold-atom optical lattices, Phys. Rev. A 73, 051601(R) (2006).
- (18) T. Matsubara, H. Matsuda, A Lattice Model of Liquid Helium, I, Prog. Theo. Phys. 16,569 (1956).
- (19) H. Matsuda, T. Tsuneto, Off-Diagonal Long-Range Order in Solids, Prog. Theo. Phys. 46,411 (1970).
- (20) For a treatment of polarized triplet on a dimer as a HCB, see K. K. Ng and T. K. Lee, Numerical study of magnetic field induced ordering in and related systems, Phys. Rev. B 73, 014433 (2006).
- (21) C. M. Varma, Missing valence states, diamagnetic insulators, and superconductors, Phys. Rev. Lett. 61, 2713 (1988).
- (22) A. Taraphder, H. R. Krishnamurthy, Rahul Pandit, and T. V. Ramakrishnan, Negative-U extended Hubbard model for doped barium bismuthates, Phys. Rev. B 52, 1368 (1995).
- (23) S. Datta, S. Yarlagadda, Supersolidity for hard-core-bosons coupled to optical phonons, Sol. State Comm. 150,2040 (2010).
- (24) G. D. Mahan, “Many Particle Physics” (Plenum Press, NY, 1990).
- (25) S. Reja, S. Yarlagadda, P. B. Littlewood, Phase diagram of the one-dimensional Hubbard-Holstein model at quarter filling, Phys. Rev. B 84,085127 (2011).
- (26) R. Pankaj and S. Yarlagadda, Study of cooperative breathing-mode in molecular chains, Phys. Rev. B 86, 035453 (2012).
- (27) R. T. Scalettar, G. G. Batrouni, A. P. Kampf, G. T. Zimanyi, Simultaneous diagonal and off-diagonal order in the Bose-Hubbard Hamiltonian, Phys. Rev. B 51,8467 (1995).
- (28) O. F. Syljuåsen, A. W. Sandvik, Quantum Monte Carlo with directed loops, Phys. Rev. E 66, 046701 (2002).
- (29) O. F. Syljuåsen, Directed loop updates for quantum lattice models, Phys. Rev. E 67, 046701 (2003).
- (30) G. G. Bartrouni, R. T. Scalettar, G. T. Zimanyi, A. P. Kampf, Supersolids in the Bose-Hubbard Hamiltonian, Phys. Rev. Lett. 74,2527 (1995).
- (31) A. W. Sandvik, Computational Studies of Quantum Spin Systems, AIP Conf. Proc. 1297, 135 (2010).
- (32) G. Schmid, S. Todo, M. Troyer, A. Dorneich, Finite-temperature Phase Diagram of Hard-Core Bosons in Two Dimensions, Phys. Rev. Lett. 88,167208 (2002).
- (33) Y. Murakami, P. Werner, N. Tsuji, and H. Aoki, Supersolid Phase Accompanied by a Quantum Critical Point in the Intermediate Coupling Regime of the Holstein Model, Phys. Rev. Lett. 113, 266404 (2014); S. Reja , S. Yarlagadda, and P. B. Littlewood, Correlated singlet phase in the one-dimensional Hubbard-Holstein model, Phys. Rev. B 86, 045110 (2012); S. Kumar and J. van den Brink, Charge ordering and magnetism in quarter-filled Hubbard-Holstein model, Phys. Rev. B 78, 155123 (2008).
- (34) Zi Cai, U. Schollwöck, and L. Pollet, Identifying a Bath-Induced Bose Liquid in Interacting Spin-Boson Models, Phys. Rev. Lett. 113, 260403 (2014).

## Acknowledgements

The authors thank Pinaki Sengupta and Keola Wierschem for useful discussions on the SSE technique and its implementation for our HCB system. We also thank Amrita Ghosh for helping with the calculations and for cross-checking.

Supplementary Material for

“Coexistence of long-range orders in a Bose-Holstein model”

Satyaki Kar and Sudhakar Yarlagadda

Comparison between ED and SSE results

We have benchmarked our QMC calculations by comparing our calculated values with those obtained by exact-diagonalization (ED) methods. We consider both the XXZ model and its simple extension, namely, the model. We find that the energy, magnetization, structure factor and SF density of our SSE calculations match quite well with those from the ED results. The comparisons of the calculated and for the model are shown in Fig. S1(a)-(b).

Our SSE results also compare well with various world-line Monte Carlo resultsscalettar1 (); scalettar2 (); scal3 () and also the SSE QMC resultschen (); wessel () available in the literature.

## References

- (35) E. Einarsson, H. J. Schulz, Direct calculation of the spin stiffness in the Heisenberg antiferromagnet, Phys. Rev. B 51,6151(1995).
- (36) K. Runge, Numerical study of the onset of superfluidity in the two-dimensional, disordered, hard-core bosons, Phys. Rev. B 45,13136(1992).
- (37) G. G. Batrouni and R. T. Scalettar, Phase Separation in Supersolids, Phys. Rev. Lett. 84, 1599 (2000).
- (38) F. Hebert, G. G. Batrouni, R. T. Scalettar, G. Schmid, M. Troyer, and A. Dorneich, Quantum phase transitions in the two-dimensional hardcore boson model, Phys. Rev. B 65, 014513 (2001).
- (39) R. T. Scalettar, G. G. Batrouni, A. P. Kampf, G. T. Zimanyi, Simultaneous diagonal and off-diagonal order in the Bose-Hubbard Hamiltonian, Phys. Rev. B 51,8467 (1995).
- (40) Y.-C. Chen, R. G. Melko, S. Wessel, Y.-J. Kao, Supersolidity from defect condensation in the extended boson Hubbard model, Phys. Rev. B 77,014524 (2008).
- (41) K.-K. Ng, Y.-C. Chen, Supersolid phases in the bosonic extended Hubbard model, Phys. Rev. B 77,052506 (2008).