{'fraction_non_alphanumeric': 0.04425638614637204, 'fraction_numerical': 0.014590189625450556, 'mean_word_length': 4.4820446735395185, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Finite temperature phases and excitations of bosons in a square lattice: A cluster mean field study

We study the finite temperature phases and collective excitations of hardcore as well as softcore bosons on a square lattice with nearest and next nearest neighbour interactions. We focus on the formation of various types of supersolid phases and their stability under thermal fluctuations. Interplay between the on-site, nearest and next nearest neighbour interactions lead to various density ordering and structural transitions, which are charted out. Thermodynamic properties and phase diagrams are obtained by extending the cluster mean field theory at finite temperatures which includes quantum effects systematically. We investigate the melting process of supersolid phase to normal fluid, which occurs in two steps due to the presence of two competing orders in the supersolid. The existence of a tetra-critical point at finite temperature exhibits intriguing behavior, which is analyzed for different regimes of interactions. We show how the phases and the transitions between them both at zero and finite temperatures can be identified from the excitation spectrum. Apart from the gapless sound mode, the appearance of Higgs like gapped mode in the supersolid phase is analyzed, which is a characteristic feature of such phase. Finally, we discuss the relevance of the results of the present work in the context of ongoing experiments on ultracold atomic gases and newly observed supersolid phases.


I. INTRODUCTION
Observation of superfluid to Mott insulator (MI) transition in a system of ultracold atoms in optical lattice [1] as well as the advancement of quantum engineering techniques has opened up the possibilities to realize various exotic phases of interacting quantum systems, which are otherwise difficult to observe in usual materials [2][3][4]. Supersolid is one of such phase which has been speculated since early 70's in the context of superfluid He 4 [5][6][7][8][9]. It is believed that, in the vicinity of solid phase, the density ordering can coexist in superfluid, resulting in a supersolid phase, which has not yet been confirmed in case of He 4 [10]. However, ultracold atomic system has become an ideal setup to study such exotic phase of matter. It turns out that in addition to the short range on-site repulsion, an effective longer range interaction is a key ingredient for supersolidity, which is present in dipolar condensate and ultracold Rydberg atoms [11][12][13][14]. Remarkably, the evidence of supersolid phase with coexisting density ordering and superfluidity has recently been confirmed experimentally in ultracold dipolar gas of 164 Dy and 168 Er atoms [15][16][17][18][19][20][21][22][23][24][25][26][27]. Apart from the dipolar condensate, supersolid phase has also been realized in a recent experiment by coupling a condensate with external cavity field [28][29][30][31][32][33]. Moreover the superstripe phase of condensate with spin orbit coupling also exhibits density modulation similar to the supersolid phase [34,35]. However, bosons in an optical lattice is an appropriate candidate to investigate the original idea of supersolidity arising as reminiscence of solid phase, where the defects of the solid can condense to form superfluid. Formation of supersolid phase in strongly correlated bosons has extensively been studied theoretically for various lattice models using different techniques [37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54]. In strongly interacting systems, particularly in one dimension, stabilization of supersolid phase with coexisting ordering is a delicate issue [55,56]. On the other hand, frustrated lattice such as triangular lattice favors the formation of supersolid [57][58][59][60][61][62][63][64][65]. Apart from the gapless sound mode which exists due to the presence of superfluidity, the supersolid phase has an additional Higgs like gapped excitation, that has been recently explored in the experiments [33]. Such excitations can serves as an important characteristic for the supersolid phase. Moreover, the nature of density ordering in supersolid can also be probed from the excitation spectrum.
In addition to quantum fluctuations, it is also a pertinent issue to investigate the the stability of supersolid phase in the presence of thermal fluctuations and its melting process at finite temperature [41,45,62,65]. In this context, the knowledge of finite temperature phase diagram also deserves further investigation to study the stability of supersolids, since the recent experiments have demonstrated the creation of such phases by directly quenching the temperature from the normal gas phase [20,26,27]. In this work, we mainly focus on the formation of supersolid phases in a system of hardcore as well as soft core bosons both at zero and finite temperatures. To incorporate the effect of long range interaction, we consider bosons with nearest (NN) and next nearest neighbor (NNN) interactions, which are crucial for the formation of supersolid phase. Moreover, the competition between NN and NNN interactions leads to different density ordering, which gives rise to additional structures of supersolid, particularly for soft core bosons. We investigate the different phases of bosons on a square lattice using finite temperature mean field theory and obtain the finite temperature phase diagram, which is important for understanding the melting of different types of solids and supersolids. Particularly the melting process of supersolid is fascinating, since it is not a priori clear which of the ordering between superfluidity and solidity will disappear first under thermal fluctuations. In addition, the characteristic features of different phases are also reflected in their excitation spectrum, which can be used for identification of respective phases as well the transition between them. For instance, the mode softening phenomena due to the appearance of 'roton' mode in homogeneous superfluid is indicative of possible density ordering leading to supersolid phase [66,67], which has also been studied experimentally [15-17, 24, 29]. Moreover different types of structural transition and density ordering can be identified from the momentum vector at which the mode softening occurs. Even the finite temperature transitions can be probed from the change in characteristic features of the excitation spectrum. Moreover, the nature of Higgs like gapped excitation in supersolid phase for this system of lattice bosons deserves further attention, since it has already been detected in trapped dipolar gas [36].
To capture the effects of correlation systematically in many body systems, mean field technique can be generalized to cluster mean field method which has been successfully applied to correlated bosons and spin systems on lattice at zero temperature [49,53,63,[68][69][70][71]. To obtain the finite temperature properties of the hardcore as well as soft core bosons with NN and NNN interaction, we extend the Cluster mean field theory at finite temperatures to capture beyond MF effects. We study how the interplay between quantum and thermal fluctuations affect the melting and ordering of various density waves and supersolid phases.
The rest of the paper is organized as follows. In Sec.II, we introduce the Hamiltonian of extended Bose Hubbard model and discuss the mean field method at zero and finite temperature. Next, the zero and finite temperature phase diagrams under the mean field approximation as well as the low lying collective excitations of different phases are discussed in Sec.III. In Sec.IV, we discuss the cluster mean field technique and investigate the different phases within this approach at finite temperatures. Finally, we conclude and summarize our results in Sec.V.

II. MODEL AND METHOD
Within tight binding approximation, bosons with a long range interaction in an optical lattice can be described by an Extended Bose Hubbard Model (EBHM), which is given by, where the different parts can be written as, The on-site part of the HamiltonianĤ given in Eq.(1) is described byĤ 0 in Eq.(2a) with the on-site interaction strength U and chemical potential µ controlling the filling of bosons. The second partĤ hop corresponds to the nearest neighbour (denoted by i, j ) hopping of the bosons with hopping amplitude t. Finally, the inter-particle interaction can be described byĤ int , consisting of nearest neighbor (NN) and next nearest neighbor (NNN) interaction (denoted by i, j ) with strength V 1 and V 2 respectively. To begin with, we review the simple mean field (MF) analysis of EBHM, which can capture the essential physics described by this model. Such interacting model can be decoupled into single site mean field HamiltonianĤ i MF described where φ i = NN j â j denotes the superfluid (SF) order parameter of i th site with j being its NN sites, whereas χ i = NN j n j and χ i = NNN j n j denote the mean field density at NN and NNN sites of i th site, respectively. In presence of only the NN interaction for square lattice, it is sufficient to consider only two site mean fields, where i ∈ A, B with A, B being the two sub-lattices shown in Fig.1(a). However, inclusion of NNN interaction induces a repulsion along the diagonal of the square lattice (see Fig.1(b-d)), which requires consideration of mean field at all the four sites for a complete description of different phases. At zero temperature, the ground state of the above Hamiltonian can be equivalently described by a Gutzwiller variational wavefunction |Ψ = i |ψ i , neglecting the inter-site correlations [72,73], where |ψ i denotes the wavefunction at each site i, Here f (n) i represents the complex variational amplitudes with the constraint n |f (n) i | 2 = 1. Within the time dependent variational theory, the dynamics can also be captured from the evolution of time dependent Gutzwiller amplitude f j | 2 , with j(j ) being summed over NN (NNN) sites of the i th site. The parameter λ i is the Lagrange multiplier at each site, associated with the constraint n |f (n) i | 2 = 1, which can be determined from the steady state configurationsf  i , the Gutzwiller amplitudes can be written as, By linearizing the above EOM in terms of the fluctuations, which can be solved in momentum space by re-writing δf n i (t) = e ı( k. ri−ω( k)t) δf (n) ( k) with r i representing the position of the i th site, we can obtain the low lying excitation spectrum ω( k) [43,51,74]. This method can also be extended to finite temperatures, where a temperature dependence is introduced in the mean fields. By diagonalizing the mean field Hamiltonian, the partition function at each site can be calculated as, At a given temperature T = 1/β, the mean field averages are computed self consistently from â i = Tr(âρ i ) and n i = Tr(nρ i ), whereρ i = e −βĤ i MF /Z i is the equilibrium density matrix. With the same spirit as the Gutzwiller method at zero temperature, we can consider the full density matrix in product form, which can be written asρ = iρ i . Also, each density matrixρ i for i th site satisfies the Liouville's equation, Within the mean field approximation, we can find the density matrix at equilibriumρ eq i fromρ eq i = 0 ⇒ [Ĥ i MF ,ρ eq i ] = 0. By imposing a small amplitude fluctuation δρ i (t) aroundρ eq i , we can perform the linear stability analysis and obtain the low lying collective modes at finite temperature [74]. In the next section, we demonstrate how such simple mean field analysis can reveal the various phases, as well the transitions among them for hardcore and softcore bosons.

III. FINITE TEMPERATURE PHASES UNDER MEAN FIELD DESCRIPTION
In this section, we elaborate the different phases of hardcore and softcore bosons in the presence of both nearest (NN) and next nearest neighbor (NNN) interaction in a square lattice and study their collective excitations within the MF approach. First, we present the MF phase diagrams charting out the various phases of this model, both at zero and finite temperatures. From the MF analysis, we obtain a Landau-Ginzburg (LG) theory of competing orders, capturing the transitions between different phases at finite temperature, as well as describing the the coexistence of superfluidity and solidity leading to supersolid (SS) phase. Finally, we discuss the low lying collective excitations of different phases at zero and finite temperatures, revealing their characteristic features as well as the structural change with increasing NNN interaction.

A. Hardcore bosons
The hardcore bosons can be realized in the limit of U → ∞, which imposes the truncation of local Hilbert space upto single occupancy of bosons, n i ∈ {0, 1}. Due to the presence of NNN interaction V 2 , we begin analyzing this system by considering a four sub-lattice structure of the unit cell of a square lattice within the MF approach. At zero temperature, the model can be described by the Gutzwiller variational wave- |1 , due to the hardcore nature of the bosons. The ground state wavefunction is obtained from minimization of energy with respect to Gutzwiller variational parameters and the different phases are identified by computing the SF order parameterφ as well as the density imbalance δ 1l = | n 1 − n l | for l ∈ {2, 3, 4} within the square lattice (cf. Fig.1). For non-zero values of δ 1l with small hopping strength t, different kinds of solid phases are formed, which can be characterised by filling n 0 . The zero temperature phase diagrams in the µ-t plane are depicted in Fig.2(a,b), which summarize all possible phases of two different regimes characterized by the NNN interaction V 2 .
First, we consider the regime 2V 2 /V 1 < 1, where the NN interaction dominates over that of NNN. As a result, around the region of half filling, a 'checkerboard' solid or density wave (DW) phase (see Fig.2(a)) is formed with vanishing SF order for small hopping amplitude t. When the hopping strength t crosses a critical value, a non vanishing SF order parameter coexist with the density ordering (δ 12 = 0), leading to the formation of checkerboard supersolid (SS CB ) phase. On the other hand when 2V 2 /V 1 > 1, the NNN interaction hinders the checkerboard ordering and a striped density ordering (δ 13 = 0) is formed in the solid and supersolid (SS STR ) phases. In addition, there are other kinds of solid phases with filling n 0 = 1/4 and 3/4, which appears below and above the half-filled region, respectively. Normally, these phases form within a narrow region of the phase diagram, however it increases with increasing the NNN strength V 2 . Subsequently, the solid phases DW-1/4 and DW-3/4 melt with increasing the hopping strength, giving rise to a new SS phase (SS3) with three sub-lattice structure (δ 12 = δ 13 = 0, see Fig.1(c)), dominated by hole or particle depending on the chemical potential µ. Further increasing the hopping strength leads to the appearance of the homogeneous SF phase with vanishing density order. Note that, the SS3 phase undergoes a first order transition to SF, whereas SS CB and SS STR undergo a second order transition to SF. We also point out the appearance of homogeneous insulating phases, namely, Mott-0 (n 0 = 0) and Mott-1 or MI (n 0 = 1), when µ is far away from the half-filled region.
Next, we discuss the various phases at finite temperature within the mean field density matrix approach, as outlined in Sec.II. For comparison with the zero temperature phase diagram, we chart out different phases in Fig.2(c,d) at finite temperatures in the µ-t plane for both the regimes of NNN interaction. Interestingly, for 2V 2 /V 1 > 1, apart from the various SS phases obtained at zero temperature, a new supersolid phase SS4 with four sub-lattice structure (δ 12 = δ 13 = δ 14 = 0, see Fig.1(d)) is formed at finite temperature in a very narrow region between the SS STR and SS3 phases, as depicted in Fig.2(d). Due to the number fluctuations at finite temperature, the insulating DW and MI phases can have non integer site occupancy with finite compressibility. As a result, for small hopping amplitudes, the different insulating phases merge together at finite T . For sufficiently large temperature, various phases shown in Fig.2(c,d), finally melt to a homogeneous normal fluid (NF) phase with vanishing SF order (see Fig.2(e,f)). We also illustrate the phase diagram in the T -t plane close to the half-filled region, summarizing the various phases for different NNN interaction strength in Fig.2(e,f). With increasing the temperature T , the DW phase melts to a normal fluid (NF) phase as the density imbalance vanishes continuously at the critical temperature T DW−NF c , leading to a second order phase transition. On the other hand, there is a continuous transition from the SF to NF phase, as the SF order vanishes at the critical temperature T SF−NF c . In the Tt plane, within the region between the DW and SF phases, the supersolid (SS) phase appears with coexisting SF order and density imbalance. Interestingly, the critical lines characterized by T DW−NF c and T SF−NF c , as well as the SS phase boundaries merge at the 'tetra-critical' point (TCP) denoted by (T * c , t * c ), which is surrounded by all four phases. Unlike the other phases, the melting of the SS phase occurs in two steps due to the coexistence of two competing orders. In the case of hardcore bosons, the superfluidity is destroyed first and the density imbalance vanishes with further increasing the temperature, resulting in the formation of the homogeneous NF. However, due to the competition between both types of orders in the SS phase, it is not straightforward to determine which type of order will be destroyed first. It is important to note that, below and close to the critical hopping strength t * c , determined by TCP, the SF can transform into SS by increasing the temperature. With increasing NNN interaction strength V 2 , the SF-SS phase boundary approaches to become vertical in the T -t plane, indicating a weak temperature dependence.

Landau-Ginzburg theory
To address the issue of competing orders, we formulate a Landau-Ginzburg (LG) theory close to the TCP, following the MF prescription. At finite temperature, under MF approximation, the density matrix of a unit cell can be written as, ρ = lρ l , where l ∈ {1, 2, 3, 4} (cf. Fig.1). Sticking close to the half-filled region, it is enough to consider a 2 sub-lattice structure with sub-lattice A(B), and the density matrix at each sub-latticeρ A(B) can be written as, Here, δ (∆) denotes the imbalance in density (SF order) between the two sub-lattices A and B. By using the density matrix in Eq.(10), we compute the average mean field (MF) free energȳ F MF as a function of m,φ, δ, and ∆ at temperature T , where S l = − ± λ ± l lnλ ± l is the von-Neumann entropy with λ ± l being the eigenvalues ofρ l . In the SS phase, due to the coexistence ofφ and δ, the SF order imbalance ∆ between the two sub-lattices also acquires a non-zero value. For the SS phase to be an equilibrium phase, we impose the steady state MF ,ρ A(B) ] = 0 in Eq. (8). As a result, we obtain ∆ in terms of δ andφ for checkerboard and striped ordering respectively as following, where, m is obtained by minimizing the free energy in Eq. (10). Hence, the free energyF MF close to TCP can be expanded in terms of δ andφ (up to fourth order) in the LG form as follows, where the coefficients are dependent on T , t and m. It should be noted that, the last term with coefficient d contributes due to the competition between the coexisting SF orderφ and density imbalance δ. By minimizing F LG with respect toφ * and δ respectively, we obtain the following linear equations, To capture the transition from DW to NF phase near the TCP, we consider the solution,φ = 0 and δ 2 = −b 2 /2b 4 . Following this, the DW-NF phase boundary is obtained from the condition b 2 = 0. Similarly, by considering the solution, δ = 0 and |φ| 2 = −a 2 /2a 4 , we obtain the SF-NF phase boundary from the condition a 2 = 0. On the other hand, to obtain the SS-DW and SS-SF phase boundaries, we consider the coexistence of ordersφ, δ = 0 and solve the linear equations in Eq.(13)(a,b), which yield the following solution, For stability of the SS phase (and consequently for existence of the tetra-critical point), it is important to have D > 0 [8]. In the inset of Fig.2(e,f), a magnified version of the immediate surrounding of TCP (pointed by a square) is illustrated, and the critical lines obtained by LG theory are shown as dotted lines. It is evident that, close to the TCP, the critical boundaries obtained from the LG theory agree well with the results obtained from numerical minimization (marked by symbols in the inset of Fig.2(e,f)).

B. Softcore bosons
For softcore bosons, the constraint on local Hilbert space is removed due to the presence of finite on-site interaction strength U , which allows multiple occupancies of bosons at each lattice site. In general, the local Hilbert space is infinitedimensional, however for the purpose of numerical analysis, we have truncated the maximum occupancy of bosons to a finite value N depending on the filling factor (density) as well on temperature. Following the similar MF procedure as discussed for the case of hardcore bosons, we obtain the different phases of softcore bosons at zero and finite temperatures.
At zero temperature, the ground state phases are obtained by minimizing the Gutzwiller variational wavefunction as discussed in Sec.II. It is to be noted that, for softcore bosons, the competition between finite on-site repulsion and NN interaction is enough to stabilize the SS phases and the presence of NNN interaction can give rise to additional structures of density ordering. First, we discuss the various phases of softcore bosons in absence of NNN interaction (V 2 = 0). In this case, the checkerboard density ordering is favoured due to NN repulsion, however the density distribution at each lattice site depends crucially on the relative strength of NN and on-site repulsion U/V 1 [43]. For sufficiently large on-site interaction U/V 1 > 4, the formation of DW phases with filling n 0 + 1/2 (where n 0 = 0, 1, . . .) occurs for reasonably small hopping strength and temperature. In addition, homogeneous Mott phases with integer filling n 0 are also formed in between two successive DW phases. On the contrary, for U/V 1 < 4, the softcore bosons prefer to stack on top of one another and also preserve the checkerboard ordering, as a consequence of which, the Mott phases disappear and are replaced by stacked checkerboard solids with filling (n 0 + 1)/2. For large hopping strength, the translational symmetry is restored, leading to the formation of homogeneous SF phase. Unlike the hardcore bosons, a SS phase is formed in between the DW and SF phase with coexisting density and SF ordering. By lowering the relative on-site interaction strength, the SS phase occupies a larger region in the phase diagram. With increasing temperature, all the phases finally melt to a homogeneous NF phase, similar to the hardcore bosons. Next, we consider the effect of NNN interaction (V 2 = 0), which leads to the additional structures in density ordering. For sufficiently large U similar to the hardcore bosons, there can be two types of density ordering, namely, checkerboard DW for 2V 2 /V 1 < 1 and striped DW for 2V 2 /V 1 > 1.
In presence of NNN interaction, the effect of on-site interaction is also relevant for density ordering. When U < min(2V 1 , 4V 2 ), the bosons accumulate on a single site, leading to the formation of new stacked DW phases with density ordering (n 0 , 0, 0, 0) in a unit cell. Note that, the configuration (n 1 , n 2 , n 3 , n 4 ) in general denotes density n l at l th site of a unit cell in the clockwise order with l = 1, 2, 3, 4 (cf. Fig.1). Interestingly, in an intermediate regime of on-site interaction, the competition between the relative strength of onsite and NNN interactions with respect to the NN interaction (U/V 1 and V 2 /V 1 respectively) leads to the formation of various types of density ordered phases, which is summarized in Fig.3 for two cases: 2V 2 /V 1 < 1 and 2V 2 /V 1 > 1. In this regime of on-site interaction, for both the cases, the homogeneous Mott phases remain absent in contrast to the case when U is strong, and instead, different types of DW phases appear, exhibiting stacking of bosons with underlying checkerboard or striped ordering. At zero temperature, for 2V 2 /V 1 < 1 and small hopping strength, such insulating phases with coexistence of checkerboard ordering as well as stacking of bosons appear for 4V 2 < U < 4(V 1 − V 2 ) (see Fig.3(a)), hence the homogeneous MI phase with unit filling is replaced by a DW phase with ordering (2, 0, 2, 0) (see Fig.4(b1,b3)). While the coexistence of stripe ordering and stacking of bosons occurs Fig.3(b)). In this case, the stacked striped DW phase with ordering (2, 2, 0, 0) appears instead of homogeneous MI phase with filling unity (see Fig.4(c1,c3)).
For U < min(2V 1 , 4V 2 ) and intermediate range of hopping strength t, a SS phase with three sub-lattice structure (SS3) is formed within a large region of phase diagram surrounding the stacked DW phases (see Fig.4(a1)). With increasing t, the SS3 phase undergoes a first order transition to the SF phase. For U > min(2V 1 , 4V 2 ) in both the regimes, 2V 2 /V 1 < 1 as well as 2V 2 /V 1 > 1, a checkerboard and striped SS phase (SS CB and SS STR ) also appear in addition to the SS3 phase, respectively (see Fig.4(b1,b3,c3)). However, as shown in Fig.4(c1), only for 2V 2 /V 1 > 1, an additional SS phase (SS4) with four sub-lattice structure can also be observed along with the SS3 and SS STR phases, when 2V 1 < U < 4V 2 . Interestingly at finite temperature, such SS4 phase also appears for U > 4V 2 within a very small region of phase diagram (see Fig.3(c4)). In general, at finite temperatures, the behavior of the different phases remains more or less similar. However, the DW phases exhibit fractional filling, and as a result, the different DW phases merge with each other for sufficiently small hopping amplitudes, which is illustrated in the bottom row of Fig.4.
We also obtain the finite temperature phase diagram in the T -t plane for fixed values of chemical potential µ around the half-filled region and for different regimes of the interaction strengths, as discussed above. For large on-site interaction strength, the phase diagram depicted in Fig.5(b) resembles to that of the hardcore bosons (see Fig.2(e)), where the different phases meet at the tetra-critical point (T * c , t * c ). On the other hand, for smaller values of U (U < 4V 2 ), the checkerboard SS phase (SS CB ) is replaced by the SS3 phase in the T -t plane and the tetra-critical point splits into a small region at finite temperature containing SS CB for 2V 2 /V 1 < 1, as depicted in Fig.5(a). Moreover, the melting of the stacked DW phase to the NF occurs via a structural transition to checkerboard solid at finite temperature. Unlike the hardcore bosons, for intermediate range of U , the two step melting of SS to NF phase can occur either through DW or SF phase, depending on the value of hopping strength t, as evident from Fig.5(c,d). The properties of the various phases such as, the structural ordering, superfluidity, as well as supersolidity, can also be reflected in the collective excitations, which we discuss in the next subsection.

C. Collective excitations
In this subsection, we investigate the low lying collective excitations of different phases of hardcore and softcore bosons at zero as well as finite temperatures following the method outlined in Sec.II.
The different features of the excitation spectrum reveals the characteristic properties of various phases, which can be used for their identification as well as transition between them. First, we investigate the excitations of the different phases at zero temperature. To understand the structural ordering, we analyze the particle hole excitations of various insulating phases at the atomic limit (t ∼ 0). The DW phase can form when the particle excitations of the vacuum state becomes unstable at µ > 0. For DW-1/4 phase, when multiple occupancy at single site is hindered by large U , two types of particle excitations are possible on the vacant sites, particularly on the diagonally opposite site and on its neighbouring site depending on the strength of V 2 with energies p3 = −µ + 4V 2 and p2 = −µ + 2V 1 , respectively. Here, pl denotes particle excitation at the l th site of a unit cell with l = 1, 2, 3, 4 (cf. Fig.1). It is evident that, for strong NN interaction strength FIG. 6: Surface plots of excitation spectrum ω( k) for hardcore bosons in the (kx/π, ky/π) plane at zero temperature for (a1,a2) V2 = 0.4 and (b1,b2) V2 = 0.6. The arrow marks in (a1,b1) for SF phases depict the roton mode softening at (1, 1) and (1, 0) in (kx/π, ky/π) plane with energy gap ∆ω(π, π) and ∆ω(π, 0) respectively. The arrow marks in (a2,b2) for SS phase point the energy gap ∆ω(π/2, π/2) and ∆ω(π/2, 0) at (1/2, 1/2) and (1/2, 0) in the (kx/π, ky/π) plane respectively.
to the checkerboard density ordering with filling n 0 = 1/2. Whereas for strong NNN interaction strength 2V 2 /V 1 > 1, the instability in p2 leads to striped ordering with filling n 0 = 1/2 when µ > 2V 1 . On the contrary, for small U i.e. U < min(2V 1 , 4V 2 ), the particle excitations at the occupied site costs a minimum energy, p1 = −µ + U , favouring a stacked DW phase for µ > U . When on-site interaction is comparable with NN and NNN interaction, it is not a priori clear, whether the particle excitations will be favourable at the unoccupied or occupied site, which can lead to the formation of different types of solid structures, coexisting with stacking. For 2V 2 /V 1 < 1, in a checkerboard solid, the particle excitation at the unoccupied sites has energy p2 = −µ + 4V 1 , whereas the creation of particle at the occupied site costs an energy p3 = −µ + 4V 2 + U . When 4V 2 < U < 4(V 1 − V 2 ), particle excitation at the occupied site becomes favourable, giving rise to a stacked checkerboard solid with configuration (2, 0, 2, 0), hindering the formation of homogeneous MI phase with unit filling, which is also evident from Fig.3(a). For similar reason, in the other case with 2V 2 /V 1 > 1, the homogeneous MI phase is absent for U < 4V 2 , due to the formation of stacked striped solids with configuration (2, 2, 0, 0). The presence of hopping strength t gives rise to the dispersion of particle hole excitations, the analytical expressions of which for the hardcore bosons are given in appendix A. The linear stability analysis can also be extended to calculate the collective excitations at finite temperatures, using Eq.(8), as discussed earlier in Sec.II. The different equilibrium phases and the transitions between them both at zero as well as finite temperature (as discussed in previous subsection), can be identified by studying the characteristic features of the low lying collective excitations. First, we focus on the collective excitations of hardcore bosons, starting from the homogeneous SF phase for different regimes of interactions. Presence of the gapless sound mode (Goldstone mode) ω( k) ∼ c s | k| with sound velocity c s signifies breaking of continuous U (1) symmetry due to non-zero SF order. In addition, near the boundary of SS phase, a 'roton' like structure is formed at higher value of momentum as a precursor of translational symmetry breaking [66,67]. For 2V 2 /V 1 < 1, such mode softening due to the formation of roton occurs at (k x /π, k y /π) = (1, 1), which leads to checkerboard density ordering. While for 2V 2 /V 1 > 1, a similar kind of soft mode appears at (k x /π, k y /π) = (1, 0), due to which it undergoes a transition to striped SS phase. For both the cases, energy gap of roton minima at corresponding momentum (indicated by arrow in Fig.6(a,b)) vanishes at the SS-SF phase boundary, indicating a continuous transition. Therefore both the SS transition and its structural transformation due to increasing NNN interactions can be detected from the mode softening of the collective excitation and its corresponding momentum.
In the SS phase, the gapless sound mode near | k| ∼ 0 persists as a result of non vanishing SF order. Moreover, due to density ordering and lattice translational symmetry breaking, the effective Brillouin zone is reduced and a gap opens up at its boundary. For 2V 2 /V 1 < 1, the gap ∆ω(π/2,π/2) at (k x /π, k y /π) = (1/2, 1/2) signifies checkerboard density ordering (see Fig.6(a2)), which vanishes at the SS-SF phase boundary. Similarly, in case of 2V 2 /V 1 > 1, a gap ∆ω(π/2, 0) at (k x /π, k y /π) = (1/2, 0) (shown in Fig.6(b2)), indicates the destruction of translational symmetry along the x-axis due to the formation of stripe phase. Hence, the structural change in the underlying sub-lattice with varying V 2 at zero temperature can be captured from the distinguishing properties of the excitation spectra, which is clearly evident from Fig.6. For density ordered phases with broken translational symmetry, the excitations can be folded into the reduced Brillouin zone and the roton branch can be viewed as a higher excitation mode with a gap, which vanishes at the SS-SF boundary. This branch of excitation in SS is similar to the amplitude (Higgs) mode of homogeneous SF at commensurate filling [75][76][77][78]. However, for clarity we show the excitation spectrum of SS phase in extended Brillouin zone.
At finite temperatures, we calculate the excitation spectrum from the fluctuations of the density matrix and obtain temperature dependence of energy gaps characterizing different phases. The excitation of NF phase is shown in Fig7(d), which exhibits a gap at k = 0 due to the absence of SF order. Melting of the SS phase can be probed from the temperature variation of both the energy gaps ∆ω(0, 0) at (k x /π, k y /π) = (0, 0) and ∆ω(π/2, π/2) at (k x /π, k y /π) = (1/2, 1/2), as shown in Fig.7(a,b). With increasing temperature, the energy gap ∆ω(0, 0) opens up at critical temperature T SS−DW c , while the ∆ω(π/2, π/2) remains finite (see Fig.7(a)), indicating a continuous transition to DW phase with vanishing SF order. Finally, as illustrated in Fig.7(b,c), the gap ∆ω(π/2, π/2) at the boundary of the reduced Brillouin zone vanishes at T DW−NF c as a result of restoration of translational symmetry in NF phase.
Finally, we discuss the features of excitation spectrum of softcore bosons in the presence of finite on-site interaction strength U for 2V 2 /V 1 < 1. In case of softcore bosons, the lowest mode exhibits characteristics akin to that of hardcore bosons, as discussed previously. However for softcore bosons, there are higher energy modes due to the enlarged local Hilbert space. Apart from the gapless Goldstone mode, there also exists a gapped Higgs mode in the SF phase Fig.8(a), which has also been observed in the recent experiments [77,78]. As discussed earlier for hardcore bosons, we noticed that a similar gapped excitation also appears in the SS phase, however the Higgs mode is absent in the SF phase. Such Higgs like gapped excitation also exist for softcore bosons, which has extensively been studied for trapped dipolar gas as a signature of SS phase [36]. Here we focus on some contrasting features of this gapped mode in SS phase (ω g−SS ) from the Higgs mode (ω H ) of the SF. In case of softcore bosons, at SF-MI transition, the gap of the Higgs mode ∆ Higgs becomes vanishingly small near the tip of the MI lobe, as shown in Fig.8(b). After entering into the insulating MI phase, the lowest excitation acquires a gap ∆ MI (see Fig.8(c)). Across the tip of the MI lobe, the SF-MI transition corresponding to the particle hole symmetry can be captured from the vanishing of the Higgs mass ∆ Higgs as well as the Mott gap ∆ MI , which is depicted in Fig.8(d). Next, we consider the scenario of SF to checkerboard SS (SS CB ) transition. In the SF phase, the gapless Goldstone mode develops a roton minima in the vicinity of the SS phase, where the Higgs mode still remains as a gapped higher energy excitation (see Fig.9(a)). The roton gap (∆ roton ) vanishes at the SF-SS boundary, indicating a continuous transition. Inside the SS phase, the Brillouin zone is reduced due to the translational symmetry breaking, consequently, a low energy gapped excitation appears apart from the higher energy excitations that correspond to the Higgs mode of SF phase, which is shown in Fig.9(b). As a signature of SS-SF continuous transition, both the roton mass in the SF and the gap of the first excitation in the SS phase (∆ g−SS ) vanish at the phase boundary (see the inset of Fig.9(b)). As a consequence of the translational symmetry breaking, the low energy Higgs like gapped excitation in the SS phase appears as a reminiscence of the roton mode, whereas the Higgs mode of SF remains as a high energy gapped excitation across the SF-SS transition. Similar behavior of the excitations can also be observed in case of striped SS (SS STR ) phase for 2V 2 /V 1 > 1along with a structural change in the underlying sub-lattice.
As discussed previously for softcore bosons, when U becomes less than the NN and NNN interactions (U < min(2V 1 , 4V 2 )), the bosons accumulate on a single site to form a stacked structure in the DW and SS phase. Similar to the previous cases, deep in the SF phase, the gapless Goldstone and gapped Higgs excitations exist, as depicted in Fig.10(a,b). When approaching the SF-SS transition, the lowest mode develops a roton minima, however the roton gap does not vanish at the phase boundary, as a consequence of a first order transition. In this regime, the SS has a 3 sublattice structure with a large density on a single site of the unit cell. As shown in Fig.10(c), the low energy excitations of SS3 consist of a gapless sound mode and a gapped mode. Finally, lowering the hopping strength leads to a stacked DW phase with a configuration (n 0 , 0, 0, 0), exhibiting the gap excitations, shown in Fig.10(d).
We have discussed the low lying excitation spectrum of hardcore and softcore bosons in the presence of NN and NNN interactions, which can be used to identify the various phases and transitions between them both at zero and finite temperature.

IV. CLUSTER MEAN FIELD THEORY
In this section, we describe the method of implementing the 'cluster mean field theory (CMFT)' [49,53,63,65,[68][69][70][71]79]. As discussed in the previous section, under MF approximation, the inter-site correlations are neglected and the total Hamiltonian of the unit cell can be decoupled individually for each site. On the contrary, in CMFT, a cluster C is considered, within which all the inter-site correlations between the bosons are treated exactly, by using exact diagonalization technique. In addition, the interaction between the bosons belonging to the given cluster C and it's neighbouring clusters are treated under mean field approximation. Therefore, the total Hamiltonian of the cluster C can be written in a more generic way as following, Where,Ĥ C describes the bosons at lattice sites within the cluster C (see Fig.11), considering the effect of correlations in an exact manner. On the other hand, the mean field interaction between the edge sites of C and its neighbouring sites outside the cluster are taken into account by the second termĤ MF , which can be written as a sum of mean field Hamiltonian for the edge sites belonging to C, ∈C n j , as well asχ i = NNN j / ∈C n j , and j (j ) denotes the NN (NNN) lattice sites outside the cluster C. However, due to the presence of sublattice symmetry throughout the entire lattice space, all the clusters can be considered to be equivalent to each other and therefore, the mean field averages . can be calculated from the lattice sites of the given cluster C. At zero temperature, the mean field averages are calculated from the ground state wavefunction |ψ by using . = ψ|.|ψ . For finite temperatures phases, the averaged mean field parameters . = Tr(.ρ T ) are obtained self-consistently by computing the equilibrium density matrixρ T = e −βĤ C /Z C at temperature T , where Z C denotes the total partition function of the cluster C. Following this procedure, the phase diagrams of hardcore as well as softcore bosons are obtained in different parameter planes, both at zero and finite temperature.
By using larger clusters and performing the finite cluster scaling [49,63], the CMFT can yield more accurate results in close agreement with that obtained from QMC, which is discussed in the appendix B. To reduce the computational cost, in the present work, we consider a 2 × 2 cluster on a two dimensional square lattice and the results are presented below.

A. Phase diagrams using CMFT
The zero temperature phase diagrams in the µ-t plane for hardcore (softcore) bosons are depicted in Fig.12(a, ( Fig.13(a1-a3)) for both the regimes of V 2 . It is evident from the figures that, all the equilibrium phases which were obtained under mean field approximation, are also present in the phase diagrams obtained by CMFT. However, there are quantitative modifications in the phase boundaries due to the presence of inter-site correlations. As a consequence of such correlations, the transition from SS/SF to DW phase occurs at smaller values of hopping strength, compared to the values obtained from the MF analysis. Furthermore, the SS regions in the phase diagrams obtained under CMFT become smaller, which were overestimated under the MF approximation. We also illustrate a set of phase diagrams in the µ-t plane for hardcore (softcore) bosons at finite temperature in Fig.12(c,d) ( Fig.13(b1-b3)). Similar to the MF description, here also the thermal fluctuations permit non integer filling at each site of the insulating phases, due to which the different solid lobes merge together when the hopping strength is sufficiently small. However, a relatively smaller region of SS phases surround them for intermediate values of hopping strength. Interestingly for hardcore bosons, unlike the MF case, the finite temperature phase diagram in the µ-t plane for 2V 2 /V 1 > 1 does not exhibit the SS4 phase in between the SS3 and SS STR regions (see Fig.12(d)), in the presence of correlations. Next, we focus on the phase diagrams of hardcore and softcore bosons in the T -t plane, obtained under CMFT, for both the regimes of V 2 , as depicted in Fig.12(e,f) and Fig.14, respectively. To compare the results obtained from CMFT with the results of MF analysis, we plot these phase diagrams close to the half-filled region. It is important to note that, due to the effect of correlations, the critical temperature corresponding to the transitions among various phases become smaller, compared to that obtained under MF description. Apart from this, in the case of hardcore bosons, another important distinguishing feature is also spotted in contrast to the MF analysis, i.e. the tetra-critical point (T * c , t * c ) observed previously (see Fig.2(e,f)), splits into two tri-critical points, (T * c , t * c ) and (T * c , t * c ) under CMFT (see Fig.12(e,f)). From a phenomenological point of view of LG theory, within the MF approach, a tetra-critical point exists due to the condition D > 0 (cf. Sec.III A 1), whereas the modification in the parameters of the LG theory incorporating the quantum fluctuations can lead to D < 0, resulting in a first order transition from DW to SF phase with the vanishing of the SS phase. Consequently, the tetra-critical point (T * c , t * c ) obtained in the MF approach transforms to a tri-critical point (T * c , t * c ). Due to the additional correlations incorporated by CMFT, there is a lowering of critical temperature T SS−DW c corresponding to the melting of SS into DW phase, and the phase boundaries containing the SS phase terminate at the second tri-critical point (T * c , t * c ) with lower temperature T * c < T * c . Consequently, the two tri-critical points (t * c , T * c ) and (t * c , T * c ) are connected by a first order phase boundary between the SF and DW phases.
For softcore bosons with finite values of on-site repulsion U , in the regime of U < 4V 2 and 2V 2 /V 1 < 1, a significant qualitative difference is noted from the phase diagram in the T -t plane (see Fig.14), when compared with the MF results (see Fig.5). As evident in Fig.14(a), in the presence of inter-site correlations, the stacked DW phase with 3 sublattice structure directly undergoes a first order phase transition into NF phase with increasing temperature, instead of melting through the checkerboard DW phase (see Fig.5(a)). Similarly, for larger values of hopping strength t, the SS CB phase observed in the MF phase diagram is not attained, when temperature of the SS3 phase is increased. Instead, the absence of SS CB phase leads to a first order phase boundary DW-SF, joining the two tri-critical points (t * c , T * c ) and (t * c , T * c ). On the other hand, when U is sufficiently large, the system approaches towards the hardcore limit, as a result, in this case also, the tetra-critical point (observed in the T -t plane within MFT in Fig.5(b)) splits into two tri-critical points under CMFT, as shown in Fig.14(b). For intermediate values of U , the nature of the tetra-critical point in the T -t plane remains as it is, although the region of the SS phase between the DW and SF phase becomes smaller compared to the MFT, which is indicated in Fig.14(c,d).
Although the CMFT significantly improves the critical points and phase boundaries obtained from MF theory, it cannot however capture the critical behavior at the vicinity of the transition. For example, capturing the Berezinskii-Kosterlitz-Thouless (BKT) transition in 2D [80] along the SF-NF phase boundary at finite temperature is beyond the scope of CMFT, as it requires much larger cluster size.

B. Superfluid fraction and compressibility
To this end, we elaborate the calculation of the different physical quantities at finite temperatures using CMFT, which are relevant for SS as well as SF phases and transition between them. Apart from the SF order parameterφ, the superfluid fraction (SFF) ρ S is also an important quantity to measure the SF transport in an interacting Bose gas. The SFF estimates the response of the system under twisted boundary conditions, which can be implemented by imposing a phase variation θ along a particular direction (say x), and the twisted HamiltonianĤ(θ) is obtained by replacing the hopping matrix element t ij with t ij e iθ in Eq.(2b), where j represents the nearest neighbour sites of i along the x direction. At finite temperatures, the SFF is defined as [81,82], where F (θ) denotes the free energy associated with the twisted HamiltonianĤ(θ), which is obtained from the partition function Z C (θ) = γ e −βEγ (θ) , with E γ (θ) being the eigenvalues ofĤ(θ) and N C denotes the number of lattice sites within the cluster C. The energies E γ (θ) can be calculated perturbatively using, whereĴ andT are the current and kinetic energy operators, respectively. These operators under CMFT are defined as follows, At zero temperature, the SFF can also be obtained from the above mentioned perturbative approach involving the opera-torsĴ andT [82]. As shown in Fig.15(b), at low temperature, with increasing hopping strength t, ρ s remains zero throughout the insulating DW phase, and becomes non-zero at the critical point t DW−SS c , indicating the appearance of a SS CB phase with non vanishing ρ S . It gradually grows with further  Fig.15(c,d)).
Next we focus on compressibility κ, which is another important physical quantity besides SFF, and is given by, Apart from SFF, the compressibility κ vanishes for insulating MI and DW phases at very low temperatures due to their incompressible nature. On the other hand, with increasing temperature, the DW/MI phases with non integer filling and the homogeneous NF phase become compressible (κ = 0) without possessing any superfluidity (ρ s = 0). Moreover, for different SS and SF phase, both κ and ρ s acquire finite values, indicating their compressible nature and superfluidity, as summarized in Fig.15(c,d). We point out that, κ exhibits a 'cusp' like structure at the transition points between different phases, which is depicted in Fig.15(b-d). In addition, close to the atomic limit (t ∼ 0), we also plot the average densityn with increasing µ (see Fig.15(a)), which exhibits a step like behavior at very low temperatures, indicating the incompressible nature of the insulating DW phases, while it is smoothened out with increasing temperature, as these phases become more compressible.
The physical quantities can also be calculated by implementing CMFT in three dimensions, although it is computationally more expensive. It is expected that the fluctuations are reduced in three dimensions, as a consequence, the results obtained from CMFT consisting of a cluster of only a single unit cell are more accurate and comparable with that of QMC, as discussed in appendix B.

V. CONCLUSION
To summarize, we have investigated the phases and collective excitations of hardcore as well as softcore bosons on a square lattice at finite temperatures, emphasizing the formation of supersolid (SS) phase and its stability under thermal fluctuations. Finite temperature properties and phase diagrams are obtained by using mean field (MF) as well as cluster mean field theory (CMFT) to study the effect of quantum fluctuations systematically. We consider bosons with nearest and next nearest neighbor interactions, which are essential for the formation of SS phase for softcore and hardcore bosons respectively. Moreover, the interplay between nearest-neighbor (NN), next nearest-neighbor (NNN) and on-site interactions leads to different types of density ordering and structures in the SS phase. We chart out the various types of density ordering in solid and SS phases, which can be categorized in three main classes. For sufficiently large on-site repulsion (and in the hardcore case), the competition between NN and NNN interactions gives rise to checkerboard and stripe density ordering. On the other hand, small on-site repulsion favors stacked structure with accumulation of particles on a single site of the unit cell. However in an intermediate regime, various types of solids can form with coexistence of stacking and checkerboard (or stripe) order.
In the phase diagram, the SS phase appears in between insulating density wave (DW) and homogeneous superfluid (SF). With increasing temperature, all the phases finally melt to the homogeneous normal fluid (NF) phase. Melting of SS is particularly interesting, since it occurs in two steps due to the presence of both the superfluidity and density ordering. Two different pathways for transformation of SS to NF phase has been analyzed from the phase diagram in the T -t plane for different parameter regimes. In T -t plane of the phase diagram, the boundaries of all the phases meet at a tetra-critical point. However from the cluster mean field results, we find for certain regime of on-site interaction, the tetra-critical point splits into two tri-critical points connected by a first order line, separating the SF and DW phases. Although the overall behavior of the MF phase diagram is similar to that obtained from CMFT, MFT over estimates the critical temperature, which can be reduced systematically by increasing the cluster size.
We have also discussed how the different phases can be identified and the transitions between them both at zero as well as finite temperatures can be probed from the excitation spectrum, which can be relevant for their experimental detection. The gapless sound mode in SF and SS signifies the presence of SF order, which vanishes in the insulating phases giving rise to a gapped excitation. In the vicinity of SS and DW phases, the excitation spectrum of SF develops a 'roton' minima, which signals possible breaking of lattice translational symmetry. Moreover, the gap between the low lying excita-tion branches at the boundary of the reduced Brillouin zone is a characteristic feature of the translational symmetry broken phases (like DW and SS), which also carries information about nature of the density ordering. At finite temperatures, we obtain the excitation spectrum from the fluctuations of the equilibrium density matrix, which can effectively capture the two step melting process from SS to NF phase. Apart from the gapless sound mode, a Higgs like gapped excitation is a characteristic feature of SS phase which has similarities with the roton mode in SF, and both become gapless at the SS-SF boundary. However for softcore bosons, the Higgs mode in SF phase does not vanish at the phase boundary.
In recent experiments, the signature of supersolidity has unambiguously been observed in trapped dipolar quantum gas of 166 Er and 164 Dy atoms [15][16][17][18][19][20][21][22][23][24][25][26][27], which shows promise to achieve the supersolid phase in strongly correlated bosons in presence of an optical lattice. The supersolid phase in trapped dipolar gas turns out to be fragile and tends to disintegrate into droplet arrays [19,20,23,24]. The presence of 2D optical lattice can stabilize such supersolid phase of dipolar atoms. Moreover, tuning the on-site interaction by Feshbach resonance can lead to the different density ordering, particularly the formation of stacked SS phase. Also the ultracold Rydberg gas is another possible candidate for supersolids, due to their long range interaction [14]. Formation of different density patterns and experimental realization of spin models with Rydberg atoms in optical lattice opens up the possibility to observe such supersolid phases of hardcore and softcore bosons [83,84]. Moreover, the tunability of effective interaction between Rydberg dressed states is suitable for the formation of supersolids with different density ordering. In recent experiments, supersolid phase of dipolar gas has been created by directly quenching the temperature, exhibiting sequential formation of density order and superfluidity [20,26,27]. These experiments show the importance to study the finite temperature phases of bosons and their out-of-equilibrium dynamics. The pathways for two step formation of SS can be probed in such quench experiments by suitably choosing the parameters. It may also be possible to observe the metastable density ordered phases in such quench process. Appearance of roton excitation as a precursor of translational symmetry broken SS phase has been extensively studied in experiments [15-17, 24, 29]. The translational symmetry broken phases (such as DW and SS) as well their melting at finite temperatures can be detected from the excitation branches, as discussed before.
In conclusion, we studied finite temperature phases and excitations of hardcore as well as softcore bosons using CMFT, focusing on the formation of SS phases with different types of density ordering. Various phases and transition between them at zero as well as finite temperature can be probed from the excitation spectrum, which has importance for their experimental detection.

Zero temperature
Following the method outlined in Sec.II, here we provide the detailed calculations and derive the analytical expressions of excitation spectrum ω( k) for different phases of hardcore bosons.
The excitation branches at zero temperature can be calculated analytically from the EOM of the time dependent Gutzwiller amplitudes, given by Eq. (5). For hardcore bosons with n i ∈ {0, 1}, the EOM for f Below, we provide the analytical expressions for dispersion of different phases of hardcore bosons at zero temperature in presence of finite hopping strength t = 0.
Checkerboard DW-1/2: By considering the two sublattice density structure (n A , n B ) ≡ (1, 0) with checkerboard ordering (cf. Fig.1(a)), the steady state values of Gutzwiller amplitudes aref (0) By substituting the steady state values in Eq.(A1) and using the steady state condition, we obtain the Lagrange multipliers at each sub-lattice, given by λ A = −µ + 4V 2 and λ B = 0. After introducing the small amplitude oscillations around the steady state valuesf i (see Eq.(6)), the linearized fluctuation equations in the momentum space within the reduced Brillouin zone can be written as following, where ( k) = 2t(cos k x + cos k y ). Solving the above equations simultaneously for ω( k) yields the particle (+) and hole (−) excitation ω ± ( k) for checkerboard DW phase with filling n 0 = 1/2, which is given by, Stripe DW-1/2: We follow the similar procedure as in the checkerboard case but with stripe ordering (cf. Fig.1(b)) in the underlying sub-lattice. From the steady state conditions of the EOM in Eq.(A1), the Lagrange multipliers at A and B sub-lattice can be fixed as λ A = −µ + 2V 1 and λ B = 0, respectively. The linearized fluctuation equations in momentum space within the reduced Brillouin zone can be written as following, where 1 ( k) = 2t cos k x , and 2 ( k) = 2t cos k y . By solving the above equations simultaneously for ω(k), the particle (+) and hole (−) excitation ω ± ( k) for striped DW phase with filling n 0 = 1/2 is given by, Homogeneous SF phase: For a homogeneous SF phase, the Gutzwiller variational wavefunction at each site can be written as |ψ = f (0) |0 +f (1) |1 , where we can use the parametrization f (0) = sin (θ/2), f (1) = cos (θ/2). From the steady state conditions of the EOM, we obtain λ = −4t cos 2 (θ/2) with θ = cos −1 µ−2(V1+V2) 2(V1+V2+2t) , and subsequently solve the following linearized fluctuation equations in the momentum space within the extended Brillouin zone, ] cos 2 (θ/2), and X( k) = V 1 (cos k x + cos k y ) + V 2 [cos(k x + k y ) + cos(k x − k y )]. The excitation spectrum ω( k) is given by, The method of linear stability analysis can be extended to finite temperatures, where the stability analysis of the finite temperature density matrixρ i at i th site is performed, which satisfies the Liouville equation in Eq. (8). We consider the following density matrix at i th site, using the above density matrix, the SF order parameter is given by, φ i = Tr(ρ iâi ) and the average density can be written in terms of m i as, n i = Tr(ρ ini ) = (1 + m i )/2. By introducing a small amplitude fluctuation δρ i (t) around the equilibrium density matrixρ eq i , such thatρ i (t) =ρ eq i +δρ i (t), and plugging it into Eq.(8), we obtain the following, where δĤ MF i contains the fluctuations around the mean field terms inĤ MF i . Substituting δρ i (t) = e ı( k. ri−ω( k)t) δρ( k) in Eq.(A12) and retaining terms only upto the first order of δρ( k), we obtain the fluctuation matrix corresponding to the set of linearized equations in momentum space. The real part of the eigenvalues ω( k) of the fluctuation matrix provides the excitation branches, whereas Im[ω( k)] = 0 ensures the stability ofρ eq . Below, we provide the semi analytical expressions for low lying excitations of some of the phases of hardcore bosons at finite temperature.
Checkerboard SS and DW: The Checkerboard SS (and DW) is an inhomogeneous phase with 2 sub-lattice structure and checkerboard ordering (cf. Fig.1(a)). In order to obtain the excitation spectrum ω( k), it requires solving a set of six fluctuation equations. The low lying excitation branches ω ± (k) for checkerboard SS at finite temperature is given by the following, where A, B, C, D are given by, , for i,ī = A, B and i =ī. Note that, here the parameters φ i and m i have an implicit dependence on temperature, which can be obtained numerically in a self consistent manner. By substituting the SF order φ i = 0 in the above equations, we obtain the excitation spectrum for checkerboard DW phase, Stripe SS and DW: Similar to the Checkerboard case, the stripe SS (and DW) is an inhomogeneous phase with 2 sub-lattice structure and stripe ordering (cf. Fig.1(b)). In this case also, we solve a set of six fluctuation equations to obtain the low lying excitations ω ± ( k) for stripe SS phase at finite temperature, which is given by the following, where A , B , C , D are given by, (A17b) , for i = A, B and i =ī. Here also, φ i and m i at finite temperature can be obtained self consistently from numerics. To get the excitation spectrum for stripe DW phase, we can substitute φ i = 0 in the above equations, which gives the following, It is easy to check, by setting m A = 1 and m B = −1 for both the Checkerboard DW (in Eq.A15) and stripe DW (in Eq.A18), the expressions match with that in Eq.(A3) and Eq.(A5) respectively, for two sub-lattice density structure (n A , n B ) = (1, 0) with filling n 0 = 1/2 at zero temperature.

Appendix B: Finite cluster size scaling of CMFT and its accuracy
Here, we check the quantitative accuracy of CMFT at finite temperature by performing the finite cluster-size scaling [49,63]. We consider the scaling parameter, λ= N B N C ×z/2 (introduced in Ref. [63]), with N B and N C denoting the number of bonds and lattice sites within the cluster C, respectively. For a square lattice, the coordination number is given by z = 4. In order to extract the critical temperature more accurately, we obtain T c for different clusters with increasing size, i.e. 2 × 2, 2 × 3, 2 × 4 (see Fig.16(a) for schematics), and finally extrapolate the values at the thermodynamic limit λ → 1 using a linear fit. For checking purposes, we consider a system of hardcore bosons, only in presence of finite NN interaction with strength V 1 , while V 2 = 0. We compare the results obtained from the CMFT at finite temperature using finite cluster-size scaling, with the Quantum Monte-Carlo (QMC) results of Ref. [85]. As shown in Fig.16(b,c), the transition points acquired from CMFT, after performing the linear extrapolation up to λ → 1, differ significantly from the single site MF result and approaches towards the more accurate value obtained from QMC. As evident from above, with the larger cluster sizes like 4 × 4 and so on, the CMFT can yield results having a better agreement with the QMC values, which can improve the phase boundaries obtained from MFT. However, working with such large clusters is challenging with our current computation power. In addition, it is not always possible to capture the true critical behavior of the transitions using finite size clusters, for e.g the SF-NF transition in 2D which is BKT like due to the absence of true long range order in SF [80], however the transition point obtained from extrapolation after using the finite size cluster scaling approaches the QMC value.
The CMFT can also be extended to three dimensions, although it costs more computation power. Here we consider DW-NF and SF-NF transitions for softcore bosons in three dimensions with NN interaction only, by considering a 2 × 2 × 2 cubic cluster, and compare the critical temperatures with that of MF and QMC results of Ref. [45]. To perform CMFT for bosons with longer range interaction can require large cluster size, which is computationally more expensive. In   17(a,b), we show the variation of density imbalance and SF order with increasing temperature, which vanishes at the critical point of DW-NF and SF-NF transition, respectively. It is evident that, the critical temperatures obtained from CMFT improves the MF value and are much closer to the QMC re-sults. In three dimensions, due to the reduced effect of fluctuations, the finite temperature CMFT results can be obtained with reasonable accuracy even with a cubic cluster of single unit cell, which is in a close agreement with the QMC values.