# Finite temperature phases and excitations of bosons on a square lattice: A
cluster mean field study

Manali Malakar  Sudip Sinha  and S. Sinha  Indian Institute of Science
Education and Research-Kolkata, Mohanpur, Nadia-741246, India

###### Abstract

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
neighbor interactions, focusing on the formation of various types of
supersolid (SS) phases and their stability under thermal fluctuations. The
interplay between the on-site, nearest, and next nearest neighbor interactions
leads to various density ordering and structural transitions, which we have
plotted out. Thermodynamic properties and phase diagrams are obtained by
cluster mean field theory at finite temperatures, which includes quantum
effects systematically, and they are compared with the single-site mean field
results. We investigate the melting process of the SS phase to normal fluid
(NF), which can occur in at least two steps due to the presence of two
competing orders in the SS. A tetra-critical point exists at finite
temperature and exhibits intriguing behavior, which is analyzed for different
regimes of interactions. The phase diagrams reveal the different pathways of
the thermal transition of SSs to the NF phase, for different interaction
regimes, which can be accessible by thermal quench protocols used in recent
experiments. We show how the phases and the transitions between them can be
identified from the characteristic features of the excitation spectrum. We
analyze the appearance of a low-energy gapped mode apart from the gapless
sound mode in the SS phase, which is analogous to the gapped mode recently
studied for dipolar SS phases. 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 SS phases.

April 2023

##  1 Introduction

Observation of superfluid to Mott insulator (MI) transition in a system of
ultracold atoms in an 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 [3, 4, 2]. Supersolid is one of such
phases which has been speculated since the late 50s in the context of
superfluid ${\rm He}^{4}$ [5, 6, 7, 8, 9, 10]. It is believed that, in the
vicinity of the solid phase, the density ordering can coexist in a superfluid,
resulting in a supersolid phase, which has not yet been confirmed in the case
of ${\rm He}^{4}$ [11]. However, the ultracold atomic system has become an
ideal setup to study such exotic phases of matter. It turns out that in
addition to the short range on-site repulsion, an effective long range
interaction is a key ingredient for supersolidity, which is present in dipolar
condensate and ultracold Rydberg atoms [12, 13, 14, 15]. Remarkably, the
evidence of a supersolid phase with coexisting density ordering and
superfluidity has recently been confirmed experimentally in ultracold dipolar
gas of ${}^{164}{\rm Dy}$ and ${}^{168}{\rm Er}$ atoms [19, 20, 21, 22, 23,
26, 27, 28, 16, 17, 18, 25, 24]. Apart from the dipolar condensate, the
supersolid phase has also been realized in a recent experiment by coupling a
condensate with an external cavity field [33, 34, 32, 29, 30, 31]. Moreover,
the superstripe phase of condensate with spin orbit coupling also exhibits
density modulation similar to the supersolid phase [35, 36]. However, bosons
in an optical lattice are an appropriate candidate to investigate the original
idea of supersolidity arising as a reminiscence of the solid phase, where the
defects of the solid can condense to form a superfluid. The formation of
supersolid phase in strongly correlated bosons has extensively been studied
theoretically for various lattice models using different techniques [53, 54,
37, 38, 39, 40, 41, 42, 43, 45, 44, 46, 47, 48, 49, 50, 55, 56, 58, 59, 61,
62, 64, 63, 51, 52, 60, 57, 65]. In strongly interacting systems, particularly
in one dimension, stabilization of the supersolid phase with coexisting
ordering is a delicate issue [66, 67, 68]. On the other hand, frustrated
lattice such as triangular lattice favors the formation of supersolid [77, 69,
70, 74, 71, 75, 76, 72, 73, 78]. The collective excitations of such phases
also exhibit interesting features. Apart from the gapless sound mode which
exists as a consequence of superfluidity, the supersolid phase has an
additional gapped excitation, that has been recently explored in the
experiments [34]. Such excitations can serve as an important characteristic of
the supersolid phase. Moreover, the nature of density ordering in supersolid
can also be probed from the excitation spectrum.

In addition to the effects of quantum fluctuations, it is also a pertinent
issue to investigate the stability of supersolid phase in the presence of
thermal fluctuations and its melting process at finite temperatures [43, 46,
49, 50, 58, 61, 74, 78]. 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
[21, 27, 28]. In this work, we mainly focus on the formation of supersolid
phases in a system of hardcore as well as softcore bosons both at zero and
finite temperatures. To incorporate the effect of long range interaction, we
consider bosons with nearest neighbor (NN) and next nearest neighbor (NNN)
interactions on a square lattice, which can play a crucial role for the
formation of supersolid phases. Moreover for softcore bosons, the interplay
between on-site, NN and NNN interactions leads to various types of density
ordering, particularly the newly observed stacked solids where bosons
accumulate on a single site of the unit cell. Such competition between the
different interactions can also give rise to additional structures of
supersolid phases. We investigate the different phases of bosons using finite
temperature mean field (MF) 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 [80, 79], which has also been studied experimentally [16, 17, 18, 25,
30]. 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. In addition, the nature of
the low energy gapped excitation in supersolid phase for this system of
lattice bosons deserves further attention, since it is analogous to the gapped
mode detected in dipolar supersolids [34, 81].

To capture the effects of correlation systematically in many body systems, the
mean field technique can be generalized to cluster mean field (CMF) method
which has been successfully applied to correlated bosons and spin systems on
lattices at zero temperature [54, 58, 61, 62, 77, 82, 83, 84, 85]. In order to
investigate the finite temperature properties of the hardcore as well as
softcore bosons with NN and NNN interaction, we extend the cluster mean field
theory at finite temperatures to capture beyond MF effects. We also 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.2, we introduce the
Hamiltonian of extended Bose Hubbard model and discuss the mean field as well
as cluster mean field method at zero and finite temperatures. Next, we study
the different phases obtained within the mean field and cluster mean field
theory at finite temperatures, as well as investigate the corresponding low-
lying collective excitations for hardcore and softcore bosons in Sec.3 and
Sec.4, respectively. In Sec.5, we summarize the stability of the various
supersolid phases obtained in the present model. Finally, we conclude and
discuss the results in Sec.6.

##  2 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,

$\displaystyle\hat{\mathcal{H}}=\hat{\mathcal{H}}_{0}+\hat{\mathcal{H}}_{\rm
hop}+\hat{\mathcal{H}}_{\rm int}$ (1)

where the different parts can be written as,

$\displaystyle\hat{\mathcal{H}}_{0}=\frac{U}{2}\sum_{i}\hat{n}_{i}(\hat{n}_{i}-1)-\mu\sum_{i}\hat{n}_{i}$
(2a) $\displaystyle\hat{\mathcal{H}}_{\rm hop}=-t\sum_{\langle
i,j\rangle}(\hat{a}^{\dagger}_{i}\hat{a}_{j}+\hat{a}^{\dagger}_{j}\hat{a}_{i})$
(2b) $\displaystyle\hat{\mathcal{H}}_{\rm int}=V_{1}\sum_{\langle
i,j\rangle}\hat{n}_{i}\hat{n}_{j}+V_{2}\sum_{\langle\langle
i,j\rangle\rangle}\hat{n}_{i}\hat{n}_{j}$ (2c)

The on-site part of the Hamiltonian $\hat{\mathcal{H}}$ given in Eq.(1) is
described by $\hat{\mathcal{H}}_{0}$ in Eq.(2a), with the on-site interaction
strength $U$ and chemical potential $\mu$ controlling the filling of bosons.
The second part $\hat{\mathcal{H}}_{\rm hop}$ corresponds to the nearest
neighbor (denoted by $\langle i,j\rangle$) hopping of the bosons with hopping
amplitude $t$. Finally, the inter-particle interaction can be described by
$\hat{\mathcal{H}}_{\rm int}$, consisting of nearest neighbor (NN) and next
nearest neighbor (NNN) interaction (denoted by $\langle\langle
i,j\rangle\rangle$) with strength $V_{1}$ and $V_{2}$ respectively.

###  2.1 Mean field theory

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
$\hat{\mathcal{H}}^{i}_{\rm MF}$ described by,

$\displaystyle\hat{\mathcal{H}}^{i}_{\rm MF}$ $\displaystyle=$
$\displaystyle-t(\hat{a}^{\dagger}_{i}\phi_{i}+\hat{a}_{i}\phi^{*}_{i})+V_{1}\hat{n}_{i}\chi_{i}+V_{2}\hat{n}_{i}\overline{\chi}_{i}$
(3) $\displaystyle+\frac{U}{2}\hat{n}_{i}(\hat{n}_{i}-1)-\mu\hat{n}_{i}$

where $\phi_{i}=\sum^{\rm NN}_{j}\langle\hat{a}_{j}\rangle$ denotes the
superfluid (SF) order parameter of $i^{\rm th}$ site with $j$ being its NN
sites, whereas $\chi_{i}=\sum^{\rm NN}_{j}\langle\hat{n}_{j}\rangle$ and
$\overline{\chi}_{i}=\sum^{\rm
NNN}_{j^{\prime}}\langle\hat{n}_{j^{\prime}}\rangle$ denote the mean field
density at NN and NNN sites of the $i^{\rm th}$ site, respectively. In
presence of only the NN interaction, for a square lattice, it is sufficient to
consider the mean field with only 2 sub-lattice structure, where $i\in A,B$
with $A,B$ being the two sub-lattices shown in Fig.1(a). However, the
inclusion of NNN interaction induces a repulsion along the diagonal of the
square lattice, which requires consideration of the mean field with a more
general 4 sub-lattice structure for a complete description of different phases
(see Fig.1(b-d)). At zero temperature, the ground state of the above
Hamiltonian can be equivalently described by a Gutzwiller variational
wavefunction $\ket{\Psi}=\prod_{i}\ket{\psi_{i}}$, neglecting the inter-site
correlations [86, 87], where $\ket{\psi_{i}}$ denotes the wavefunction at each
site $i$,

$\displaystyle\ket{\psi_{i}}=\sum_{n}f^{(n)}_{i}\ket{n_{i}}$ (4)

Figure 1: Schematics of the possible underlying sub-lattice structures in a
unit cell for different equilibrium phases of bosons in a square lattice. For
the 2 sub-lattice structure in (a,b), the sub-lattice $A\,(B)$ is denoted by
black (white) circles.

Here $f^{(n)}_{i}$ represents the complex variational amplitudes with the
constraint $\sum_{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^{(n)}_{i}(t)$ as follows,

$\displaystyle\imath\dot{f}^{(n)}_{i}\\!=\\!-t\left(\phi^{*}_{i}\sqrt{(n+1)}f^{(n+1)}_{i}+\phi_{i}\sqrt{n}f^{(n-1)}_{i}\right)$
$\displaystyle+\\!f^{(n)}_{i}\left[\frac{U}{2}n\left(n-1\right)-\left(\mu-
V_{1}\chi_{i}-V_{2}\bar{\chi}_{i}\right)n-\lambda_{i}\right]$ (5)

where, $\phi_{i}=\sum_{j,n}\\!\sqrt{(n+1)}f^{(n+1)}_{j}\\!f^{*(n)}_{j}$,
$\chi_{i}=\sum_{j,n}n|f^{(n)}_{j}|^{2}$, and
$\bar{\chi}_{i}=\sum_{j^{\prime},n}n|f^{(n)}_{j^{\prime}}|^{2}$, with
$j(j^{\prime})$ being summed over NN (NNN) sites of the $i^{\rm th}$ site. The
parameter $\lambda_{i}$ is the Lagrange multiplier at each site, associated
with the constraint $\sum_{n}|f^{(n)}_{i}|^{2}=1$, which can be determined
from the steady state amplitudes $\bar{f}^{(n)}_{i}$, describing the
equilibrium phases. In presence of small time dependent fluctuations $\delta
f^{(n)}_{i}(t)$ around the steady state values $\bar{f}^{(n)}_{i}$, the
Gutzwiller amplitudes can be written as,

$\displaystyle f_{i}^{(n)}(t)=\bar{f}^{(n)}_{i}+\delta f^{(n)}_{i}(t)$ (6)

By linearizing the above equations of motion (EOM) in terms of the
fluctuations, which can be solved in momentum space by re-writing $\delta
f^{(n)}_{i}(t)=e^{\imath(\vec{k}.\vec{r}_{i}-\omega(\vec{k})t)}\delta
f^{(n)}(\vec{k})$ with $\vec{r}_{i}$ representing the position of the $i^{\rm
th}$ site, we can obtain the low-lying excitation spectrum $\omega(\vec{k})$
[45, 57]. This method can also be extended to finite temperatures, where a
temperature dependence is introduced in the mean fields [78]. From the mean
field Hamiltonian in Eq.(3), the partition function at each site can be
calculated as,

$\displaystyle Z_{i}={\rm Tr}\left(e^{-\beta\hat{\mathcal{H}}^{i}_{\rm
MF}}\right)$ (7)

At a given temperature $T=1/\beta$, the mean field averages are computed self
consistently from $\langle\hat{a}_{i}\rangle={\rm Tr}(\hat{a}\hat{\rho}_{i})$
and $\langle\hat{n_{i}}\rangle={\rm Tr}(\hat{n}\hat{\rho}_{i})$, where
$\hat{\rho}_{i}=e^{-\beta\hat{\mathcal{H}}^{i}_{\rm MF}}/Z_{i}$ is the
equilibrium density matrix at the $i^{\rm th}$ site. 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
$\hat{\rho}=\prod_{i}\hat{\rho}_{i}$ [88]. Also, each density matrix
$\hat{\rho}_{i}$ satisfies the Liouville’s equation,

$\displaystyle\dot{\hat{\rho}}_{i}=-\imath[\hat{\mathcal{H}}^{i}_{\rm
MF},\hat{\rho}_{i}]$ (8)

Within the mean field approximation, we can find the equilibrium density
matrix $\hat{\rho}^{\rm eq}_{i}$ from $\dot{\hat{\rho}}^{\rm
eq}_{i}=0\Rightarrow[\hat{\mathcal{H}}^{i}_{\rm MF},\hat{\rho}^{\rm
eq}_{i}]=0$. By imposing a small amplitude fluctuation
$\delta\hat{\rho}_{i}(t)$ around $\hat{\rho}^{\rm eq}_{i}$, we can perform the
linear stability analysis and obtain the low-lying collective modes at finite
temperatures. Next, we consider the effect of correlations among the different
sites by implementing the ‘cluster mean field (CMF) theory’ and explore the
different phases of hardcore and softcore bosons beyond the MF analysis both
at zero as well as finite temperatures.

###  2.2 Cluster mean field theory

Figure 2: Schematics of a $2\times 2$ cluster for extended Bose Hubbard Model
(EBHM), which is denoted by the solid blue square. The dotted red (dashed
black) lines denote the NN (NNN) mean field interactions with the lattice
sites which are not within the cluster $\mathcal{C}$.

In this section, we describe the method of implementing the ‘cluster mean
field (CMF) theory’ [54, 58, 61, 62, 77, 78, 82, 83, 84, 85, 89]. 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 CMF approach, a cluster
$\mathcal{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
edge sites of the given cluster $\mathcal{C}$ and its neighboring clusters are
treated under mean field approximation. Therefore, the total Hamiltonian of
the cluster $\mathcal{C}$ can be written in a more generic way as following,

$\displaystyle\hat{\mathcal{H}}=\hat{\mathcal{H}}_{\mathcal{C}}+\hat{\mathcal{H}}_{\rm{MF}}$
(9)

where, $\hat{\mathcal{H}}_{\mathcal{C}}$ describes the bosons at lattice sites
within the cluster ${\mathcal{C}}$ (see Fig.2), considering the effect of
correlations in an exact manner. On the other hand, the MF interaction between
the edge sites of $\mathcal{C}$ and its neighboring sites outside the cluster
are taken into account by the second term $\hat{\mathcal{H}}_{\rm{MF}}$, which
can be written as a sum of the mean field Hamiltonian for the edge sites
belonging to $\mathcal{C}$,

$\displaystyle\hat{\mathcal{H}}_{\rm{MF}}$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\sum_{i\,\in\,{\rm
edge\,sites}}\\!\\!\\!\\!\hat{\mathcal{H}}^{i}_{\rm{MF}}$ (10)
$\displaystyle\hat{\mathcal{H}}^{i}_{\rm MF}$ $\displaystyle=$
$\displaystyle-t(\hat{a}^{\dagger}_{i}\phi_{i}+\hat{a}_{i}\phi^{*}_{i})+V_{1}\hat{n}_{i}\chi_{i}+V_{2}\hat{n}_{i}\overline{\chi}_{i}$
(11)

where $\phi_{i}=\sum^{\rm NN}_{j\notin\mathcal{C}}\langle\hat{a}_{j}\rangle$,
$\chi_{i}=\sum^{\rm NN}_{j\notin\mathcal{C}}\langle\hat{n}_{j}\rangle$, as
well as $\bar{\chi}_{i}=\sum^{\rm
NNN}_{j^{\prime}\notin\mathcal{C}}\langle\hat{n}_{j^{\prime}}\rangle$, and $j$
($j^{\prime}$) denote the NN (NNN) lattice sites outside the cluster
$\mathcal{C}$. However, due to the presence of sub-lattice symmetry throughout
the entire lattice space, all the clusters can be considered to be equivalent
to each other and therefore, the MF averages $\langle.\rangle$ can be
calculated from the lattice sites of the given cluster $\mathcal{C}$. At zero
temperature, the mean field averages are obtained from the ground state
wavefunction $|\psi\rangle$ by using
$\langle.\rangle=\langle\psi|.|\psi\rangle$. For finite temperature phases,
the thermal averages $\langle.\rangle$ = $\mathrm{Tr}(.\>\hat{\rho}_{T})$ are
obtained self-consistently by computing the equilibrium density matrix
$\hat{\rho}_{T}=e^{-\beta\hat{\mathcal{H}}_{\mathcal{C}}}/\mathcal{Z}_{\mathcal{C}}$
at temperature $T$, where $\mathcal{Z}_{\mathcal{C}}$ denotes the total
partition function of the cluster $\mathcal{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 temperatures, which are
discussed in the subsequent sections.

By using larger clusters and performing the finite cluster-size scaling [54,
77], the CMF theory can yield more accurate results in close agreement with
that obtained from Quantum Monte-Carlo (QMC) simulations, which is discussed
in the A. To reduce the computational cost, in the present work, we provide
the results by considering a $2\times 2$ cluster in the unit cell of a two
dimensional square lattice, which will be discussed in the next section.

##  3 Hardcore Bosons

In this section, we revisit the formation of different phases of hardcore
bosons in presence of both the nearest neighbor (NN) and next nearest neighbor
(NNN) interaction on a square lattice. Here, our main motivation is to study
the finite temperature phases beyond the MF analysis by systematically
incorporating the correlations using the cluster mean field theory. We also
study the collective excitations within the MF approach, which can be used to
identify the various phases and the transition between them at finite
temperatures.

###  3.1 Finite temperature Phases

To begin with, we chart out various phases of hardcore bosons both at zero and
finite temperatures, within the MF and CMF theory. The main motivation of the
present analysis is to study the different phases, particularly the formation
of supersolid (SS) phases, due to the competition between NN and NNN
interactions at finite temperatures. Hardcore bosons can be realized in the
limit of large on-site interaction ($U\rightarrow\infty$), which imposes the
constraint on local Hilbert space up to single occupancy at each lattice site,
$n_{i}\in\\{0,1\\}$. We compare the phase diagrams obtained within the single
site MF theory with that of CMF theory for cluster size
$N_{\mathcal{C}}=2\times 2$ in the unit cell of a square lattice (see Fig.2).
At zero temperature, within the single site MF approach, the ground state can
be described by the Gutzwiller variational wavefunction,

$\displaystyle\ket{\Psi}=\prod_{i}\left(f^{(0)}_{i}\ket{0}+f^{(1)}_{i}\ket{1}\right),$
(12)

where the variational parameters $f^{(n)}_{i}$ are determined from
minimization of energy, subjected to the constraint
$|f^{(0)}_{i}|^{2}+|f^{(1)}_{i}|^{2}=1$. To explore the rich variety of phases
in presence of both $V_{1}$ and $V_{2}$, we consider an unit cell with 4 sub-
lattice structure (cf. Fig.1), where the different density ordering can be
characterized from the density imbalance
$\delta_{1l}=|\langle\hat{n}_{1}\rangle-\langle\hat{n_{l}}\rangle|$ (for
$l\in\\{2,3,4\\}$) within the unit cell. The presence of superfluidity is
determined by the non vanishing average SF order parameter
$\bar{\phi}=\sum^{4}_{l=1}\phi_{l}/4$.

As a result of the competition between NN ($V_{1}$) and NNN ($V_{2}$)
interactions, two different types of density ordering can be obtained around
the half-filled region ($n_{0}=1/2$): checkerboard (CB) ordering for
$2V_{2}\\!<\\!V_{1}$ and striped (STR) ordering for $2V_{2}\\!>\\!V_{1}$ [39].
In the regime, where NN interaction dominates over that of NNN i.e.
$2V_{2}\\!<\\!V_{1}$, it is expected that with the increasing hopping strength
$t$, a ‘checkerboard’ supersolid can form around the solid phase. However,
from the previous Quantum Monte-Carlo (QMC) studies, it has been shown that
such supersolid phase at zero temperature is thermodynamically unstable and
undergoes a phase separation [40]. Although the checkerboard SS phase can be
spuriously obtained from the MF theory, by incorporating the correlations
within the CMF theory, its region of formation diminishes with increasing the
cluster size [54]. However, such checkerboard SS phase can be stabilised in
presence of a realistic long range interaction [53, 54] or by using a super-
lattice structure [59].

On the other hand, for $2V_{2}\\!>\\!V_{1}$, it is known that a stable stripe
supersolid phase SSSTR can emerge with coexisting SF order ($\bar{\phi}\neq
0$) and stripe density ordering ($\delta_{13}\neq 0$), surrounding the solid
phase [39, 40, 54]. Due to this reason, we mainly focus on the finite
temperature phases in this regime in order to study the stability of the SS
phases and their melting under thermal fluctuations. First, we plot out the
different phases at zero temperature in this regime obtained from the MF and
the CMF method, as depicted in the phase diagram of Fig.3(a). Apart from the
striped solid with filling $n_{0}=1/2$, other solids with filling $n_{0}=1/4$
(DW-1/4) and $n_{0}=3/4$ (DW-3/4) also appear in the hole and particle
dominated regions, respectively, which give rise to a SS phase (SS3) with 3
sub-lattice structure ($\delta_{12}\neq\delta_{13}\neq 0$, see Fig.1(c)), as
the hopping strength $t$ increases. Within the CMF theory with $2\times 2$
cluster, the region of SS phase in the phase diagram decreases due to the
incorporation of correlations. Furthermore, increasing the hopping strength
leads to the appearance of the homogeneous SF phase.

Next, we discuss the various phases at finite temperatures using the density
matrix formalism within the MF and CMF approaches, as outlined in Sec.2. For
comparison with the zero temperature phase diagram, we obtain the different
phases at finite temperatures in the $\mu/V_{2}$-$t/V_{2}$ plane, which is
shown in Fig.3(b). At finite temperatures, the insulating phases (DW and Mott
insulators) can have non integer site occupancies with vanishing SF order,
however, thermally enhanced number fluctuations give rise to finite
compressibility (see Sec.4.1 for more details in case of softcore bosons). In
spite of becoming compressible at finite temperatures, these phases remain
insulating as the transport coefficients characterized by the superfluid
fraction (SFF) vanish (see for example Fig.7(b) and Sec.4.1 for more details
on SFF). Although the finite temperature DW phases have non-integer filling at
different lattice sites, their density ordering is still retained at low
temperatures and the different insulating lobes merge together for small
hopping strengths. Interestingly, due to such non-integer occupancies at
different lattice sites, an insulating DW3 phase with a 3 sub-lattice
structure ($\delta_{12}\neq\delta_{13}\neq 0$) can appear in the particle and
hole dominated regions, above and below the striped ordered DWSTR phase,
respectively (see Fig.3(b)). With further increasing the temperature, these
phases undergo a transition to a homogeneous normal fluid (NF) phase. On the
other hand, the insulating Mott phases at finite temperature have a
homogeneous structure with non-integer filling, which exhibit a smooth
crossover to the NF phase [89, 90]. At a sufficiently high temperature, all
the phases depicted in Fig.3(b) finally melt to the NF phase with vanishing SF
order (see Fig.3(c)). As seen from Fig.3(a,b), although the MF analysis
captures all the equilibrium phases, the effect of correlations within the CMF
theory can modify the phase boundaries with more quantitative accuracy (see
also A). Consequently, the supersolid regions become narrower compared to that
obtained from the single site MF theory.

Figure 3: Mean field (MF) and cluster mean field (CMF) phase diagrams for
hardcore bosons in different parameter planes, in the regime
$2V_{2}\\!>\\!V_{1}$ ($V_{1}/V_{2}=0.44$). Various phases are depicted in the
$\mu/V_{2}$-$T/V_{2}$ plane: (a) at zero and (b) at finite temperature
$T/V_{2}=0.05$. Note that, DW3 and SS3 have a 3 sub-lattice structure (cf
Fig.1(c)). (c) Phase diagram in the $T/V_{2}$-$t/V_{2}$ plane for fixed value
of chemical potential $\mu/V_{2}=1.8$. The tetra-critical point
$(T^{*}_{c}/V_{2},t^{*}_{c}/V_{2})$ within MF theory is denoted by a solid
arrow, and the tri-critical points
$(T^{*^{\prime}}_{c}\\!\\!/V_{2},t^{*^{\prime}}_{c}\\!/V_{2})$ as well as
$(T^{*^{\prime\prime}}_{c}\\!\\!\\!/V_{2},t^{*^{\prime\prime}}_{c}\\!\\!/V_{2})$
within CMF theory are denoted by dashed arrows. (d) Finite temperature phases
in the $T/V_{2}$-$V_{1}/V_{2}$ plane for $t/V_{2}=0.25$ and $\mu/V_{2}=0.8$.
In these and the rest of the phase diagrams, the thick blue solid/dashed lines
represent the phase boundaries within the MF theory, whereas the red
open/closed squares represent the same within CMF theory. The thick blue solid
lines (red open squares) denote second order transition, and thick blue dashed
lines (red closed squares) denote first order transition within MF (CMF)
theory. Note that, the red thin solid/dashed lines act only as guide to eye
for CMF phase boundaries.

We also obtain the phase diagram in the $T/V_{2}$-$t/V_{2}$ plane close to the
half-filled region for a fixed value of chemical potential $\mu$, using both
the MF and CMF approach. As seen from Fig.3(c), with increasing temperature,
the $\rm DW_{\rm STR}$ phase melts to the NF phase as the density imbalance
vanishes continuously at the critical temperature $T^{\rm DW-NF}_{c}/V_{2}$,
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 $\bar{\phi}$
vanishes at the critical temperature $T^{\rm SF-NF}_{c}/V_{2}$. In the
$T/V_{2}$-$t/V_{2}$ plane, the SSSTR phase appears within the region between
the DWSTR and SF phases. Unlike the other phases, the melting of the SS phase
occurs in two steps due to the coexistence of two competing orders, where the
superfluidity is destroyed first and the density imbalance vanishes with
further increasing the temperature, resulting in the formation of the
homogeneous NF. Similar two step thermal transition of SS phase has also been
observed in previous QMC studies [46, 50]. Interestingly, under MF
approximation, the critical lines characterized by $T^{\rm DW-NF}_{c}/V_{2}$
and $T^{\rm SF-NF}_{c}/V_{2}$, as well as the SS phase boundaries merge at the
‘tetra-critical’ point (TCP) denoted by ($T^{*}_{c}/V_{2},t^{*}_{c}/V_{2}$),
which is surrounded by all four phases (see Fig.3(c)). The presence of such
TCP can be explained phenomenologically from Landau Ginzburg theory of
competing orders [9, 49]. As mentioned earlier, the region of SS phase shrinks
due to the effect of correlations within the CMF theory, which leads to a
distinct feature of splitting of the TCP $(T_{c}^{*}/V_{2},t_{c}^{*}/V_{2})$
into two tri-critical points,
$(T_{c}^{*^{\prime}}\\!\\!/V_{2},t_{c}^{*^{\prime}}\\!\\!/V_{2})$ and
$(T_{c}^{*^{\prime\prime}}\\!\\!/V_{2},t_{c}^{*^{\prime\prime}}\\!/V_{2})$,
depicted in Fig.3(c). In the CMF phase diagram, there is a lowering of
critical temperature $T_{c}^{\rm SS-DW}/V_{2}$ corresponding to the melting of
${\rm SS}_{\rm STR}$ to the ${\rm DW}_{\rm STR}$ phase, and the phase
boundaries of the supersolid phase terminate at the second tri-critical point
($T^{*^{\prime\prime}}_{c}\\!\\!/V_{2},t^{*^{\prime\prime}}_{c}\\!\\!/V_{2}$)
with lower temperature $T^{*^{\prime\prime}}_{c}\\!<\\!T_{c}^{*^{\prime}}$.
Consequently, the two tri-critical points
$(T^{*^{\prime}}_{c}\\!\\!/V_{2},t^{*^{\prime}}_{c}\\!/V_{2})$ and
$(T^{*^{\prime\prime}}_{c}\\!\\!/V_{2},t^{*^{\prime\prime}}_{c}/V_{2})$ are
connected by a first order line between the SF and ${\rm DW}_{\rm STR}$
phases.

Next, we explore the finite temperature phases and the structural transitions
arising due to the competing NN interaction strength $V_{1}$, for a fixed
chemical potential $\mu$ around the quarter-filled region $(n_{0}=1/4)$, which
is depicted in the phase diagram in $T/V_{2}$-$V_{1}/V_{2}$ plane for moderate
values of hopping strength (see Fig.3(d)). With increasing the interaction
strength $V_{1}$, a structural transition in both the SS and DW phases is
observed. The striped supersolid SSSTR undergoes a first order transition to
the SS3 phase with a 3 sub-lattice structure, and for higher temperatures, the
DWSTR phase transforms to the DW3 phase. Such structural transition has also
been observed in QMC studies [48, 50]. On the contrary, within the single site
MF theory, such discontinuous transition is smeared out at finite temperatures
by the appearance of DW4 and SS4 phases (with a 4 sub-lattice structure) in a
narrow region between the DWSTR/SSSTR and DW3/SS3 phases, as shown in the
inset of Fig.3(d).

###  3.2 Collective excitations

In this subsection, we study the low-lying collective excitations of different
phases of hardcore bosons at zero and finite temperatures, to identify the
characteristic features of different phases which are charted out in the
previous subsection. We calculate the collective excitations at zero
temperature from the small amplitude fluctuations of time-dependent Gutzwiller
wavefunction, which has also been extended to obtain the finite temperature
excitations using the fluctuations of density matrix in Eq.(8), as outlined in
Sec.2.

Figure 4: Excitation spectrum of hardcore bosons obtained from MF theory at
zero temperature for (a) superfluid (SF) and (b) striped supersolid (SSSTR)
phase, in the regime $2V_{2}\\!>\\!V_{1}$ $(V_{1}/V_{2}=0.44)$. The Brillouin
zone in (b) is folded along $k_{x}$, whereas the inset shows excitation within
the extended scheme. Various energy gaps are denoted by double headed arrows.
The transitions from (c) supersolid (SSSTR) to DWSTR and (d) DWSTR to normal
fluid (NF) phases with increasing temperature $T/V_{2}$, are illustrated
through the variation of excitation gaps (order parameters): $\Delta(0,0)$
($\bar{\phi}$) and $\Delta(\pi/2,0)$ ($\delta_{13}$), respectively. In the
inset of (d), the behavior of the energy gap $\Delta(\pi/2,0)$ is zoomed close
to the critical temperature $T^{\rm DW-NF}_{c}/V_{2}$ corresponding to the
continuous DW-NF transition. Parameters chosen: $\mu/V_{2}=1.8$, and for
(c,d): $t=0.844t^{\rm SS-SF}_{c}$. The transition point in (a,b) is $t^{\rm
SS-SF}_{c}/V_{2}\sim 0.474$.

To this end, we mainly focus on the NNN dominated regime
($2V_{2}\\!>\\!V_{1}$) to study the various characteristic features of the
collective excitations associated with the formation of SS phase at zero and
finite temperatures. First, we calculate the excitation spectrum of the
homogeneous SF phase in the vicinity of stripe supersolid (SSSTR), around the
half-filled region. The presence of the gapless sound mode (Goldstone mode)
$\omega(\vec{k})\sim c_{s}|\vec{k}|$ with sound velocity $c_{s}$ signifies
breaking of continuous $U(1)$ symmetry in SF phase. In addition, as shown in
Fig.4(a), near the boundary of SSSTR phase, a ‘roton’ mode arises at momentum
vector $(k_{x},k_{y})=(\pi,0)$ and $(0,\pi)$ with gap $\Delta_{\rm roton}$, as
a precursor of translational symmetry breaking [79, 80], which vanishes at the
SF-SS boundary, indicating a second order transition to SSSTR phase.
Subsequently, the striped density ordering can occur along the $x$ or
$y$-axis. For the formation of stripes along the $y$-axis, a gap opens up
along $k_{x}=\pi/2$ and the Brillouin zone is reduced within this region of
$k_{x}$ (See Fig.4(b)), with minimum gap at $(k_{x},k_{y})=(\pi/2,0)$. Such
gap $\Delta(\pi/2,0)$ characterizes the stripe ordering in both the SS and DW
phases. For clarity, we have also shown the excitation along $k_{x}$ in the
extended Brillouin zone in the inset of Fig.4(b). However, the SS phase
exhibits a gapless sound mode due to the presence of non vanishing SF order
parameter, whereas a gap $\Delta(0,0)$ opens up at $\vec{k}=0$ in the
insulating DW phase. As a consequence of folding the Brillouin zone, an
additional excitation branch appears in the SSSTR phase with a gap
$\Delta_{\rm g-SS}$ at $\vec{k}=0$, which vanishes at the boundary of the
continuous SS-SF transition.

The above features of excitations not only distinguish the different phases,
but the variation of different energy gaps with interaction strength as well
with temperature, can also be utilized to identify the zero and finite
temperature transitions. As shown in Fig.4(c), with increasing the
temperature, the opening of the gap $\Delta(0,0)$ of the lowest excitation
branch signals the SS-DW transition, that correctly captures the critical
temperature $T^{\rm SS-DW}_{c}\\!/V_{2}$ obtained from the MF theory, above
which the average SF order parameter $\bar{\phi}$ vanishes. The striped
density ordering of the SS/DW phase is characterized from the non vanishing
energy gap $\Delta\left(\pi/2,0\right)$, which can also be used to capture the
transition from the homogeneous phases to the phases with broken translational
symmetry. The melting of DWSTR phase to a homogeneous NF phase with increasing
temperature (see Fig.3(c)) is also signalled as the energy gap
$\Delta\left(\pi/2,0\right)$ closes at the critical temperature $T^{\rm DW-
NF}_{c}\\!/V_{2}$, which is shown in Fig.4(d). In addition, from the variation
of the energy gap $\Delta\left(\pi/2,0\right)$ and density imbalance
$\delta_{13}$ with temperature, we observe an appearance of a kink at $T^{\rm
SS-DW}_{c}\\!/V_{2}$, corresponding to the SS-DW transition. Such features of
the excitation spectrum obtained from the time dependent density matrix method
are consistent with the mean field results, which can also be relevant for
experimental detection of finite temperature phases. We also provide the
analytical and semi-analytical expressions of excitation spectrum for
different phases of hardcore bosons at zero and finite temperatures,
respectively, in B. In the next section, we study the finite temperature
phases and collective modes of the softcore bosons in presence of NN and NNN
interactions.

##  4 Softcore Bosons

In this section, we study the additional density ordered phases arising in
case of softcore bosons at finite temperatures. 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 our numerical studies, we truncate the maximum occupancy
of bosons at each site to a finite value $N$ depending on the filling factor
(density) and temperature. Following the similar procedure as discussed for
the case of hardcore bosons, we explain different phases of softcore bosons
obtained under mean-field approximation and compare the results with CMF
theory, which includes the effects of correlation. Finally, we identify
different features of these phases as well as transitions between them by
analysing the characteristic properties of the excitation spectrum, at zero
and finite temperatures.

Figure 5: The mean field phase diagram in $\mu/V_{1}$-$U/V_{1}$ plane
illustrating the different insulating (DW) phases close to the atomic limit
($t/V_{1}\sim 0$), for (a) $V_{2}/V_{1}=0.4$ and (b) $V_{2}/V_{1}=0.6$. The
vertical dotted lines marked by the solid arrows represent the critical
boundaries of the different regimes of interaction strengths. (c-h) Schematics
of the different possible density wave (DW) phases. (c,f) denote density
ordering with bosons occupied on a single site of the unit cell, (d,g)
represent checkerboard like density ordering, whereas (e,h) indicate that of
striped like. Note that, the configuration $(n_{1},n_{2},n_{3},n_{4})$ in
general denotes density $n_{l}$ at $l^{\rm th}$ site of the unit cell in the
anti-clockwise order with $l=1,2,3,4$.

###  4.1 Finite temperature Phases

Before studying the finite temperature behavior of softcore bosons, we discuss
various phases at zero temperature, also focusing on their density ordering,
which arises as a result of the competition between the on-site ($U$), NN
($V_{1}$), and NNN ($V_{2}$) interaction strengths. Unlike the hardcore
bosons, the NN interaction with strength $V_{1}$ is sufficient for stabilising
the SS phase for softcore bosons [37, 44], whereas the presence of NNN
interaction ($V_{2}$) can give rise to a rich variety of phases with
additional structures of density ordering.

Figure 6: MF and CMF phase diagram for softcore bosons in the
$\mu/V_{1}$-$T/V_{1}$ plane at zero (top row) and at finite temperature
$T/V_{1}=0.05$ (bottom row). Here, the supersolid phase SS4 in (c,f) has a 4
sub-lattice structure (cf. Fig.1(d)). Note that, the DW3 (stacked) phases in
(d) have a 3 sub-lattice structure (cf. Fig.3) which has a large density on a
particular lattice site, and is obtained from the melting of $(n_{0},0,0,0)$
solid phases with $n_{0}\geq 2$. Also, the DW3 phases (with 3 sub-lattice
structure) in (e,f) at higher values of $\mu$ have the same order of
occupancies as their zero temperature counterparts (in (b,c)), but with a non-
integer filling at each site. Parameters chosen for (a,b), (d,e):
$V_{2}/V_{1}=0.4$, and for (c,f): $V_{2}/V_{1}=0.6$.

To begin with, we review 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 the NN repulsion, however the density distribution within
the unit cell depends crucially on the relative strength of on-site repulsion
with respect to NN interaction i.e. $U/V_{1}$ [44, 55, 45]. 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,\ldots$), having the configuration
$(n_{0}+1,n_{0},n_{0}+1,n_{0})$, occurs for reasonably small hopping strength
and temperature. Note that, the configuration $(n_{1},n_{2},n_{3},n_{4})$ in
general denotes density $n_{l}$ at $l^{\rm th}$ site of a unit cell in the
anti-clockwise order with $l=1,2,3,4$ (cf. Fig.1). In addition, homogeneous
Mott phases with integer filling $n_{0}$ can also form 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, also preserving the checkerboard
ordering, where the nearest neighbor of the filled sites always remain empty
[45]. As a consequence, the Mott phases disappear and are replaced by stacked
checkerboard solids with filling $n_{0}/2$, having the configuration
$(n_{0},0,n_{0},0)$. 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.

Next, we incorporate the effect of NNN interaction ($V_{2}\neq 0$), in the
presence of which, a richer variety of different density ordered phases
appear, as charted out in the $\mu/V_{1}$-$U/V_{1}$ plane in Fig.5, that can
be discussed in three different regimes of interaction strengths.

Case-I: For $U\\!\gtrsim\\!{\rm min}(2V_{1},4V_{2})$ and small hopping
strength, two types of density ordered phases with filling $n_{0}=1/2$ can
appear, namely: checkerboard solid for $2V_{2}\\!<\\!V_{1}$
($U\\!\gtrsim\\!4V_{2}$), and striped solid for $2V_{2}\\!>\\!V_{1}$
($U\\!\gtrsim\\!2V_{1}$). Moreover, for sufficiently large $U$, along with the
above mentioned solids, the homogeneous Mott insulator (MI) phase with unit
filling ($n_{0}=1$) can appear in the regime $U\\!\gtrsim\\!4(V_{1}-V_{2})$
for $2V_{2}\\!<\\!V_{1}$ and $U\\!\gtrsim\\!4V_{2}$ for $2V_{2}>V_{1}$,
similar to the hardcore bosons.

Case-II: On the contrary, when $U\\!\lesssim\\!{\rm min}(2V_{1},4V_{2})$, the
bosons accumulate on a single site, leading to the formation of new stacked
solid phases with density ordering $(n_{0},0,0,0)$ in a unit cell.

Case-III: Interestingly, in an intermediate regime of on-site interaction
strength $U$ (in between case-I and II), its competition with relative
strength $V_{2}/V_{1}$, leads to the formation of various other types of
density ordered phases, as shown in Fig.5 for both $2V_{2}\\!<\\!V_{1}$ and
$2V_{2}\\!>\\!V_{1}$. In this regime, the homogeneous Mott phases remain
absent in contrast to the case-I when $U$ is sufficiently strong, 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}$ and small hopping strength, such insulating phases with
coexistence of checkerboard ordering as well as stacking of bosons appear for
$4V_{2}\\!\lesssim\\!U\\!\lesssim\\!4(V_{1}-V_{2})$ (see Fig.5(a)). Here, the
homogeneous MI phase with unit filling (appearing in case-I for sufficiently
large $U$) is replaced by a DW phase with configuration $(2,0,2,0)$, while the
coexistence of stripe ordering and stacking of bosons occurs for
$2V_{2}\\!>\\!V_{1}$ when $2V_{1}\\!\lesssim\\!U\\!\lesssim\\!4V_{2}$ (see
Fig.5(b)). In this case, the stacked striped DW phase with configuration
$(2,2,0,0)$ appears instead of the unit filled MI phase. It is important to
note that, in presence of both $V_{1}$ and $V_{2}$, such phenomena of stacking
with an underlying checkerboard ordering can also lead to the formation of new
solid phases with 3 sub-lattice structure, having a configuration of
$(n_{0}+1,0,n_{0},0)$ (see Fig.5(a,g) and 6(b)), which can appear in between
two successive insulating lobes with density ordering $(n_{0},0,n_{0},0)$ and
$(n_{0}+1,0,n_{0}+1,0)$, within the NN dominated regime
($2V_{2}\\!<\\!V_{1}$). On the other hand, similar behavior of stacking is
observed but with a striped ordering in the NNN dominated regime
($2V_{2}\\!>\\!V_{1}$), where the new solid phases have a configuration of
$(n_{0}+1,n_{0},0,0)$ (see Fig.5(b,h) and 6(c)), which appear in between
$(n_{0},n_{0},0,0)$ and $(n_{0}+1,n_{0}+1,0,0)$ insulating lobes.

To this end, we discuss how the above mentioned solids (DW) with different
density ordering lead to the formation of corresponding SS phases by
increasing the hopping strength $t$. The phase diagrams obtained within the MF
and CMF theory are depicted in $\mu/V_{1}$-$t/V_{1}$ plane of Fig.6, for
different regimes of on-site interaction strength $U$. For
$U\\!\lesssim\\!{\rm min}(2V_{1},4V_{2})$ (case-II), a SS phase with 3 sub-
lattice structure (SS3) is formed in a large region of the phase diagram for
intermediate range of hopping strength $t$, surrounding the stacked solid
phases (see Fig.6(a)). As $t$ increases, the SS3 phase undergoes a first order
transition to the homogeneous SF phase. In addition to the SS3 phase, for
intermediate values of $U$ (case-III), when
$4V_{1}\\!\lesssim\\!U\\!\lesssim\\!4(V_{1}-V_{2})$ for $2V_{2}\\!<\\!V_{1}$
and $2V_{1}\\!\lesssim\\!U\\!\lesssim\\!4V_{2}$ for $2V_{2}\\!>\\!V_{1}$, a
significant region of checkerboard as well as striped SS phase (${\rm
SS_{CB}}$ and ${\rm SS_{STR}}$) appear, respectively (see Fig.6(b,c)), and
both become gradually smaller with increasing $U$ (case-I). However, as shown
in Fig.6(c), within NNN dominated regime $2V_{2}\\!>\\!V_{1}$, an additional
SS phase (SS4) with a 4 sub-lattice structure can also be observed along with
the SS3 and ${\rm SS_{STR}}$ phases, around the particle dominated region for
larger values of the chemical potential $\mu$.

Next, we compare the finite temperature phases obtained from the MF and CMF
theory with the zero temperature phase diagrams for different regimes of
interaction strengths in Fig.6(d-f). We have observed that all the zero
temperature phases are retained at low temperatures and an additional
homogeneous NF phase can appear at higher temperatures. It is evident from the
phase diagrams in Fig.6 that, all the equilibrium phases obtained from the MF
theory are also present within the CMF approach, however, it leads to some
quantitative modifications. Notably, the effect of correlations can
significantly reduce the SS regions, in between the SF and DW phases, which is
overestimated in the results obtained from single-site MF approximation. In
addition, at finite temperatures, the thermal fluctuation leads to the
formation of insulating phases with non integer filling, due to which the
insulating lobes with commensurate fillings merge together when the hopping
strength is sufficiently small, as seen from Fig.6(d-f), even though the
respective density orderings within the insulating lobes are still preserved
at low temperatures. To illustrate this, we compute the average density
$\bar{n}=\sum_{i\in N_{\mathcal{C}}}\langle\hat{n}_{i}\rangle/N_{\mathcal{C}}$
and the corresponding number fluctuations $\Delta n=\sum_{i\in
N_{\mathcal{C}}}(\langle\hat{n}^{2}_{i}\rangle-\langle\hat{n}_{i}\rangle^{2})/N_{\mathcal{C}}$
within the cluster $\mathcal{C}$, with increasing $\mu$ for two different
values of temperature, close to the atomic limit ($t\sim 0$) in the regime
corresponding to Fig.6(e). As noted from Fig.7(a), the average density
$\bar{n}$ exhibits a step like behavior at low temperatures, indicating the
incompressible nature of the insulating DW phases, while it is smeared out
with increasing $T$, due to the thermally enhanced number fluctuations.
Although these phases (DW and Mott insulators) become compressible at finite
temperatures, they remain insulating as the transport coefficients quantified
from the superfluid fraction (SFF) vanish (see Fig.7(b)).

Figure 7: (a) Variation of average density $\bar{n}$ (left axis) and
corresponding number fluctuations $\Delta n$ (right axis) with increasing
$\mu$ at different temperatures for $t/V_{1}=0.05$. (b) Behavior of the
different physical quantities: density imbalance $\delta_{12}$, average SF
order $\bar{\phi}$, and superfluid fraction (SFF) $\rho_{s}$ are illustrated
with increasing hopping strength $t/V_{1}$, at a fixed temperature using CMF
theory. The vertical dashed lines denote the transition points, separating the
different phases. (c) Comparison of SF fraction $\rho_{s}$ obtained within MF
(solid line) and CMF (dotted line) at zero temperature. The vertical solid
(MF) and dotted (CMF) lines indicate the DWCB-SSCB (dark magenta lines) and
SSCB-SF (light orange lines) transition. In the inset, the different
contributions of SFF
$\rho_{s}=\rho^{\text{\tiny(1)}}_{s}-\rho^{\text{\tiny(2)}}_{s}$ (see text for
details) arising from $\rho^{\text{\tiny(1)}}_{s}$ (thick lines) and
$\rho^{\text{\tiny(2)}}_{s}$ (thin lines) are shown for MF (solid lines) and
CMF (dotted lines), respectively. Parameters chosen: $\mu/V_{1}=2.5$,
$U/V_{1}=2.0$ and $V_{2}/V_{1}=0.4$.

Next, we explore the finite temperature behavior of the superfluid fraction
(SFF) $\rho_{s}$, which serves as another important physical quantity apart
from the average SF order parameter $\bar{\phi}$, that quantifies the
superfluid transport property of an interacting Bose gas. The SFF estimates
the response of the system under twisted boundary conditions along a
particular direction (say $x$), so that the twisted Hamiltonian
$\hat{\mathcal{H}}(\theta)$ is obtained by replacing the nearest neighbor
hopping matrix element $t_{ij}$ with $t_{ij}e^{i\theta}$
($\theta\\!\ll\\!\pi$) in Eq.(2b). At finite temperatures, the SFF is defined
as [95, 96],

$\displaystyle\rho_{s}=\frac{F(\theta)-F(0)}{N_{\mathcal{C}}t\theta^{2}}\>\>;\>F(\theta)=-T\>\mathrm{ln}\mathcal{Z}_{\mathcal{C}}(\theta)$
(13)

where $F(\theta)$ denotes the free energy associated with the twisted
Hamiltonian $\hat{\mathcal{H}}(\theta)$, which is obtained from the partition
function $\mathcal{Z}_{\mathcal{C}}(\theta)$ = $\sum_{\gamma}e^{-\beta
E_{\gamma}(\theta)}$, with $E_{\gamma}(\theta)$ being the eigenvalues of
$\hat{\mathcal{H}}(\theta)$ and $N_{\mathcal{C}}$ denotes the number of
lattice sites within the cluster $\mathcal{C}$. At zero temperature, the free
energy in Eq.14 can be replaced by the ground state energy to obtain the SFF.
The energies $E_{\gamma}(\theta)$ in presence of the twist can be calculated
perturbatively using,

$\displaystyle E_{\gamma}(\theta)$ $\displaystyle=$ $\displaystyle
E_{\gamma}(0)+\theta^{2}\sum_{\gamma^{{}^{\prime}}\neq\gamma}\frac{|\langle\gamma|\hat{\mathcal{J}}|\gamma^{{}^{\prime}}\rangle|^{2}}{E_{\gamma}(0)-E_{\gamma^{{}^{\prime}}}(0)}-\frac{\theta^{2}}{2}\langle\gamma|\hat{\mathcal{T}}|\gamma\rangle$
(14)

where $\ket{\gamma}$ ($\gamma=0,1,2...$) denotes the unperturbed eigenstates
with corresponding eigenvalues $E_{\gamma}(0)$ of the Hamiltonian
$\hat{\mathcal{H}}(0)$ in absence of twist. Here, $\hat{\mathcal{J}}$ and
$\hat{\mathcal{T}}$ are the current and kinetic energy operators,
respectively. These operators under CMF theory are defined as follows [89],

$\displaystyle\hat{\mathcal{J}}\\!=\\!-\imath
t\Big{[}\sum_{i<j\in\mathcal{C}}\hat{a}^{\dagger}_{i}\hat{a}_{j}+\\!\\!\\!\\!\sum_{\langle
i,j\rangle,i\in\mathcal{C},j\notin\mathcal{C}}\hat{a}^{\dagger}_{i}\langle\hat{a}_{j}\rangle\Big{]}+\mathrm{h.c}$
(15a)
$\displaystyle\hat{\mathcal{T}}\\!=\\!-t\Big{[}\sum_{i<j\in\mathcal{C}}\hat{a}^{\dagger}_{i}\hat{a}_{j}+\\!\\!\\!\\!\sum_{\langle
i,j\rangle,i\in\mathcal{C},j\notin\mathcal{C}}\hat{a}^{\dagger}_{i}\langle\hat{a}_{j}\rangle\Big{]}+\mathrm{h.c}$
(15b)

where $j$ corresponds to the NN sites of the $i^{\rm th}$ site along a
particular direction (say $x$-axis). At zero temperature, the SFF can also be
obtained from the above mentioned perturbative approach involving the
operators $\hat{\mathcal{J}}$ and $\hat{\mathcal{T}}$. Using the second order
perturbation method, the SFF can be written as
$\rho_{s}=\rho^{\text{\tiny(1)}}_{s}-\rho^{\text{\tiny(2)}}_{s}$, where the
contributions from the kinetic energy and current operators are given by,
$\rho^{\text{\tiny(1)}}_{s}=-\langle
0|\hat{\mathcal{T}}|0\rangle/2tN_{\mathcal{C}}$ and
$\rho^{\text{\tiny(2)}}_{s}=\sum_{\gamma\neq
0}|\langle\gamma|\hat{\mathcal{J}}|0\rangle|^{2}/\left[tN_{\mathcal{C}}(E_{\gamma}(0)-E_{0}(0))\right]$,
respectively [96]. As shown in Fig.7(b,c), with increasing hopping strength
$t$, $\rho_{s}$ becomes non vanishing above the critical point $t^{\rm DW-
SS}_{c}/V_{1}$, manifesting the transition from DWCB to SSCB. Further
increasing the hopping strength, $\rho_{s}$ exhibits a sharp change,
indicating the SS-SF transition. In Fig.7(c), we also compare the SFF as well
as its different components, $\rho^{\text{\tiny(1)}}_{s}$ and
$\rho^{\text{\tiny(2)}}_{s}$, obtained from MF and CMF theory at zero
temperature, to show the effect of correlations within the CMF approach.

Figure 8: Finite temperature phases of softcore bosons in the
$T/V_{1}$-$t/V_{1}$ plane within MF and CMF theory for fixed values of
chemical potential $\mu$. In (a,b), both the MF and CMF phases are depicted
together, whereas in (c,d), the MF counterparts are shown separately in the
inset. The tetra-critical points $(T^{*}_{c}/V_{1},t^{*}_{c}/V_{1})$ in
(a,b,d) and $(\tilde{T}^{*}_{c}/V_{1},\tilde{t}^{*}_{c}/V_{1})$ in (d) are
marked by the solid arrows, whereas the tri-critical points
$(T^{*^{\prime}}_{c}\\!\\!/V_{1},t^{*^{\prime}}_{c}\\!/V_{1})$ and
$(T^{*^{\prime\prime}}_{c}\\!\\!/V_{1},t^{*^{\prime\prime}}_{c}\\!\\!/V_{1})$
in (c) are marked by the dotted arrows. Note that, the DW3 (stacked) phase in
(c) is obtained from the melting of $(2,0,0,0)$ solid phase, whereas the DW3
phase in (d) is obtained from the melting of $(2,0,1,0)$ solid phase.

The physical quantities can also be calculated by implementing CMF theory 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 a cluster consisting of only a single
unit cell are more accurate and comparable with that of QMC, as discussed in
A.

To focus on the melting of different DW and SS phases formed for various
regimes of interaction strengths, we also obtain the finite temperature phase
diagrams in the $T/V_{1}$-$t/V_{1}$ plane for fixed values of chemical
potential $\mu$ around the half-filled region. As discussed previously, for
sufficiently large on-site repulsion $U$, two types of density ordering can
occur near the half-filled region, which include checkerboard order for
$2V_{2}\\!<\\!V_{1}$ and striped order for $2V_{2}\\!>\\!V_{1}$. The melting
of such ordered phases with increasing temperature are summarized in Fig.8(a)
and 8(b), respectively. It is evident from Fig.8(a,b), there exists a tetra-
critical point (TCP) denoted by $(T_{c}^{*}/V_{1},t_{c}^{*}/V_{1})$, which is
surrounded by the different phases. Within CMF theory, the presence of
correlations not only lowers the values of
$(T_{c}^{*}/V_{1},t_{c}^{*}/V_{1})$, but also significantly shrinks the SS
regions. On the other hand, for smaller values of on-site interaction
$U\\!\lesssim\\!{\rm min}(2V_{1},4V_{2})$, in the regime of
$2V_{2}\\!<\\!V_{1}$, the checkerboard ordering in DW and SS phases is
hindered, and instead, the stacking of bosons on a single site leads to the
formation of stacked DW3 and SS3 phases, respectively. Within the single site
MF theory, unlike the previous cases, the TCP is replaced by a small region
containing the SSCB phase at finite temperatures, as depicted in the inset of
Fig.8(c). Moreover, with increasing the temperature, the melting of stacked
DW3 phase to the homogeneous NF phase occurs in two steps via a first-order
structural transition to DWCB phase. However, due to the inclusion of
correlations, in the phase diagram obtained within CMF theory (see Fig.8(c)),
such melting process occurs directly from the DW3 to NF phase via a first
order transition, even though, the ${\rm DW}_{\rm CB}$ phase can occur as a
metastable state at higher energy. Due to the competing interaction strengths,
different types of density ordered phases can exist as metastable states [47],
which can be probed through the non-equilibrium processes, typically by
quenching the temperature [63, 21, 27, 28]. It is interesting to note, the
region of SSCB phase at finite temperature obtained from the MF theory
disappears in the CMF phase diagram, resulting in the appearance of two tri-
critical points $(T_{c}^{*^{\prime}}\\!\\!/V_{1},t_{c}^{*^{\prime}}\\!/V_{1})$
and
$(T_{c}^{*^{\prime\prime}}\\!\\!/V_{1},t_{c}^{*^{\prime\prime}}\\!\\!/V_{1})$,
joined by a first order phase boundary, separating the DW3 and SF phases, as
shown in Fig.8(c). In addition, we also study the melting of DW-$(2,0,1,0)$
phase (having a 3 sub-lattice structure) with increasing temperature and
hopping in the $T/V_{1}$-$t/V_{1}$ plane in Fig.8(d), away from half filling
in the particle dominated region. It is observed that, as the temperature
increases, both the stacking and checkerboard ordering coexist
($\delta_{13}\neq\delta_{12}\neq 0$) in the DW3 phase till a certain critical
temperature, above which only the checkerboard ordering prevails with a 2 sub-
lattice structure ($\delta_{13}=0$, $\delta_{12}\neq 0$), resulting in the
formation of DWCB, that ultimately melts to the homogeneous NF phase
($\delta_{13}=\delta_{12}=0$). Unlike the previous phase diagrams around the
half-filled region corresponding to Fig.8(a,b), in this case, an additional
tetra-critical point (TCP) denoted by
($\tilde{T}^{*}_{c}/V_{1},\tilde{t}^{*}_{c}/V_{1}$) is observed, with
$\tilde{T}^{*}_{c}<T^{*}_{c}$ and $\tilde{t}^{*}_{c}<t^{*}_{c}$, where the
different DW/SS phases with 2 and 3 sub-lattice structure meet together. Here,
the two TCPs are joined by a second order phase boundary, separating the DWCB
and SSCB phases. For comparison purpose, the MF counterpart of this phase
diagram is also illustrated in the inset of Fig.8(d).

Figure 9: Two step melting of checkerboard supersolid (SSCB) with increasing
temperature via the two possible pathways in the regime corresponding to
Fig.8(a). Variation of the different physical quantities: density imbalance
$\delta_{12}$, average SF order $\bar{\phi}$, and superfluid fraction (SFF)
$\rho_{s}$ are shown with increasing temperature at a fixed hopping strength
(c) below and (d) above the tetra-critical value $t^{*}_{c}/V_{1}\simeq
0.183$, using CMF theory. The vertical dashed lines indicate the transition
points, separating the different phases. Parameters chosen: $\mu/V_{1}=2.5$,
$U/V_{1}=2.0$, and $V_{2}/V_{1}=0.4$.

Similar to the hardcore bosons, with increasing temperature, the melting of
the SS phase can occur in at least two steps due to the coexistence of
competing orders i.e. density modulation ($\delta_{1l}$) and SF order
($\bar{\phi}$). Interestingly, as seen from Fig.8(a,b), for moderate values of
$U$, both types of SS phases formed for $2V_{2}\\!<\\!V_{1}$ and
$2V_{2}\\!>\\!V_{1}$, can melt to the NF phase through two possible pathways
across the phase boundaries located at either side of the TCP
$(T_{c}^{*}/V_{1},t_{c}^{*}/V_{1})$. For $t<t^{*}_{c}$, the two step melting
of the SS phase occurs through the formation of DW phase, which finally melts
to NF phase. On the other hand, the second pathway of melting of SS can occur
when $t>t^{*}_{c}$, where the density ordering vanishes first, leading to the
formation of SF phase, which finally undergoes a transition to NF phase. The
two step melting process of the SSCB phase in the regime $2V_{2}\\!<\\!V_{1}$
is illustrated from the variation of different physical quantities in
Fig.9(a,b). It is important to note that, for very large values of $U$, as the
hardcore limit is approached, in both the regimes of $V_{2}$ mentioned above,
the SS-SF phase boundary shown in Fig.8(a,b) bends towards the left side, thus
eliminating the possibility of SS-SF-NF transition pathway (see for example
Fig.3(c)). Interestingly, for small on-site interaction strength
$\left(U\\!\lesssim\\!{\rm min}(2V_{1},4V_{2})\right)$, a three step melting
process of SS3 phase via SS3-DW3-SF-NF can occur when
$t^{*^{\prime}}_{c}<t<t^{*^{\prime\prime}}_{c}$, as evident from Fig.8(c).
Moreover, a similar three step thermal transition SS3-DW3-DWCB-NF can also
occur in the regime corresponding to Fig.8(d) when $t^{\rm DW-
SS}_{c}<t<\tilde{t}^{*}_{c}$.

As mentioned previously, although the CMF approach significantly improves the
critical points and phase boundaries obtained from the MF theory, it is
difficult to capture the critical behavior at the vicinity of the transition.
For example, capturing the Berezinskii-Kosterlitz-Thouless (BKT) transition in
2D [97] along the SF-NF phase boundary at finite temperature is beyond the
scope of CMF theory, as it requires much larger cluster size. 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.

###  4.2 Collective excitations

In this subsection, we investigate the low-lying collective excitations of
different phases of softcore bosons, following the method outlined in Sec.2.
The different features of the excitation spectrum reveal the characteristic
properties of various phases, which can be used for their identification as
well as transitions between them. We particularly investigate the excitations
of the different phases at zero temperature, the properties of which remain
mostly similar at finite temperatures, although the various energy gaps can
vary with temperature.

Figure 10: (a-c) Few low-lying excitations of softcore bosons are shown at
zero temperature along the diagonal of the Brillouin zone
$k_{x}\\!=\\!k_{y}\\!=\\!k_{d}$, for different values of hopping strength $t$
along the commensurate filling $n_{0}=1$. The gapless (Higgs) mode in SF phase
is denoted by solid (dashed) line in (a,b). The gap of Higgs mode $\Delta_{\rm
Higgs}$ in SF phase and lowest mode $\Delta_{\rm MI}$ in MI phase at $k_{d}=0$
are marked by a double-headed arrow in (a) and (c), respectively. (d)
Variation of the Higgs mass $\Delta_{\rm Higgs}$ (solid line) in the SF phase
and the Mott gap $\Delta_{\rm MI}$ (dotted line) in the MI phase with hopping
strength $t$. Parameters chosen: $\mu/V_{1}=7.25$, $U/V_{1}=4$,
$V_{2}/V_{1}=0.4$. The transition point is $t^{\rm MI-SF}_{c}/V_{1}\simeq
0.171$

To understand the structural ordering of different insulating phases, we first
analyze their particle hole excitations in the atomic limit ($t\sim 0$). The
DW phase in the presence of NN interaction can form when the particle
excitations of the vacuum state becomes unstable at $\mu\\!>\\!0$. For large
values of $U$, the DW-1/4 phase (single particle at $l=1$ site of the unit
cell) can have two possible particle excitations on the vacant sites,
particularly on the diagonally opposite site, which costs an energy
$\epsilon_{p3}=-\mu+4V_{2}$, whereas in the neighboring site, its energy is
given by $\epsilon_{p2}=-\mu+2V_{1}$. Here, $\epsilon_{pl}$ denotes the
particle excitation energy at the $l^{\rm th}$ site of a unit cell with
$l=1,2,3,4$ (cf. Fig.1). It is evident that, for strong NN interaction
strength ($2V_{2}\\!<\\!V_{1}$), $\epsilon_{p3}$ becomes unstable first when
$\mu\\!>\\!4V_{2}$, leading to the checkerboard density ordering with filling
$n_{0}=1/2$. On the other hand for $2V_{2}\\!>\\!V_{1}$, the instability in
$\epsilon_{p2}$ occurs first when $\mu\\!>\\!2V_{1}$, leading to a striped
solid with filling $n_{0}=1/2$. In contrast, for small $U$ i.e. $U\\!<\\!{\rm
min}(2V_{1},4V_{2})$, the particle excitations at the occupied site costs a
minimum energy, $\epsilon_{p1}=-\mu+U$, favouring a stacked DW phase for
$\mu\\!>\\!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 solids of different density ordering coexisting with stacking.
For $2V_{2}\\!<\\!V_{1}$, in a checkerboard solid, the particle excitation at
the unoccupied sites has energy $\epsilon_{p2}=-\mu+4V_{1}$, whereas the
creation of particle at the occupied site costs an energy
$\epsilon_{p3}=-\mu+4V_{2}+U$. When $4V_{2}\\!<\\!U\\!<\\!4(V_{1}-V_{2})$, the
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.5(a). For similar reason, in the other case with
$2V_{2}\\!>\\!V_{1}$, the homogeneous MI phase is absent for
$2V_{1}\\!<\\!U\\!<\\!4V_{2}$, due to the formation of stacked striped solids
with configuration $(2,2,0,0)$ (see Fig.5(b)). Such simple analysis of the
particle hole excitations of the insulating phases in the atomic limit can
give an idea of the possible density orderings of the softcore bosons.
Moreover, these density orderings can give rise to structures in the SS phase,
which finally disappear when the hopping strength or temperature becomes
sufficiently large, leading to the formation of homogeneous phases. With
finite hopping strength, such particle hole excitations acquire a dispersion,
which we investigate for different regimes of interactions.

Figure 11: Few low-lying excitations of softcore bosons are shown at zero
temperature along the diagonal of the Brillouin zone
$k_{x}\\!=\\!k_{y}\\!=\\!k_{d}$, for different values of hopping strength
$t/V_{1}$ near the SSCB phase within the half filled region. The gapless
(Higgs) mode is denoted by solid (dashed) line in (a). The Brillouin zone is
folded along the diagonal $k_{d}$ in (b) and the white headed arrow points
towards the energy gap ($\Delta_{\rm g-SS}$) of the lowest gapped mode (blue
dashed dotted line) at $k_{d}=0$. The inset in (b) shows the variation of the
roton gap $\Delta_{\rm roton}$ (solid line) in the SF phase and the gap
$\Delta_{\rm g-SS}$ (dotted line) in the SSCB phase, with increasing $t$.
Parameters chosen: $\mu/V_{1}=2.6$, $U/V_{1}=4$, and $V_{2}/V_{1}=0.4$. The
transition point is $t^{\rm SS-SF}_{c}/V_{1}\simeq 0.25$. Figure 12: Few low-
lying excitations of softcore bosons are shown at zero temperature along the
diagonal of the Brillouin zone $k_{x}=k_{y}=k_{d}$, for different values of
hopping strength $t/V_{1}$ near the SS3 phase around half filled region. The
gapless (Higgs) mode is denoted by solid (dashed) line in (a,b). The Brillouin
zone is folded along the diagonal $k_{d}$ in (c,d). Parameters chosen:
$U/V_{1}=1.3$, $V_{2}/V_{1}=0.4$, and $\mu/V_{1}=2$. The transition points are
$t^{\rm SS-SF}_{c}/V_{1}\simeq 0.273$ in (a,b) and $t^{\rm DW-
SS}_{c}/V_{1}\simeq 0.17$ in (c,d).

We mainly focus on the excitation spectrum in the regime $2V_{2}\\!<\\!V_{1}$,
where the checkerboard structure is favoured for $U\gtrsim 4V_{2}$, as
discussed previously. In the SF phase, apart from the gapless Goldstone mode,
there also exists a gapped mode which is the usual Higgs mode in the SF phase
arising due to the continuous $U(1)$ symmetry breaking [91, 92], as shown in
Fig.10(a), that has recently been observed in experiments [93, 94]. It is
important to note that, for SF-MI continuous transition, the gap of the Higgs
mode $\Delta_{\rm Higgs}$ vanishes near the tip of the MI lobe close to the
commensurate filling $n_{0}=1$, as shown in Fig.10(b), which is consistent
with the previous observations [92]. On the other hand, in the MI phase, the
particle hole excitations acquire an energy gap at $\vec{k}=0$ (see
Fig.10(c)). The variation of the lowest excitation gap in the Mott phase
($\Delta_{\rm MI}$) as well as the Higgs mass $\Delta_{\rm Higgs}$ in the SF
phase with increasing hopping strength is depicted in Fig.10(d), which
captures the second order SF-MI transition across the tip of the MI lobe.
Next, we consider the scenario of SF to checkerboard SS (SSCB) transition. As
discussed previously, the SSCB phase surrounds the checkerboard solid (DWCB)
phase, so that the density imbalance in the two sub-lattices ($\delta_{12}$)
changes continuously, and finally transforms into a homogeneous SF. It is
evident from Fig.11(a) that, in the SF phase, the gapless mode (at
$\vec{k}=0$) develops a roton minima with energy gap $\Delta_{\rm roton}$ at
$k_{x}\\!=\\!k_{y}\\!=\\!\pi$, which vanishes at the SF-SS phase boundary,
signalling a continuous transition to the translational symmetry-broken SSCB
phase. The appearance of such mode softening phenomena at the aforementioned
momentum vector $\vec{k}=(\pi,\pi)$ is indicative of the checkerboard density
ordering. In the SSCB phase, in addition to the gapless sound mode arising due
to non-vanishing SF order, there also exist two gapped modes. For the higher
energy gapped mode, the energy gap remains finite and changes continuously
across the SS-SF transition (see Fig.11), after which, this mode transforms
into the Higgs mode of the SF phase, as discussed before. On the other hand,
the low-lying gapped excitation in the SSCB phase (with energy gap
$\Delta_{\rm g-SS}$) is a reminiscent of the roton mode of the SF phase, which
emerges due to the folding of the Brillouin zone, as a gap in the excitation
opens at $k_{x}\\!=\\!k_{y}\\!=\\!\pi/2$, due to the lattice translational
symmetry breaking. It is interesting to note that, an analogous low energy
gapped excitation has also been observed in the SS phase of a dipolar gas [81,
34], which can be identified as a Higgs mode arising due to the breaking of
continuous translation symmetry. However, in the present case, no continuous
translational symmetry breaking is associated with the formation of the SS
phases. As shown in the inset of Fig.11(b), both the energy gaps: $\Delta_{\rm
g-SS}$ of the low-lying gapped mode of SSCB and $\Delta_{\rm roton}$ of the SF
phase, vanish at the critical point $t^{\rm SS-SF}_{c}\\!/V_{1}$, reflecting
the continuous SS-SF transition. Similar behavior of the excitations can also
be observed in the case of striped SS (SSSTR) to SF phase transition for
$2V_{2}\\!>\\!V_{1}$ and $U\gtrsim 2V_{1}$, which has been discussed
previously in details for the case of hardcore bosons in Sec.3.2.

As discussed before for softcore bosons, when $U\\!\lesssim\\!{\rm
min}(2V_{1},4V_{2})$, the bosons accumulate on a single site to form a stacked
structure in the DW and SS phases. In this regime, the SS phase (SS3) has a 3
sub-lattice structure with a large density on a single site of the unit cell.
Due to the 3 sub-lattice structure, the density imbalances within the unit
cell ($\delta_{12}$ and $\delta_{13}$) do not change continuously to form a
homogeneous SF phase, indicating a first order SS-SF transition. Note that,
the presence of $V_{2}$ is important for the formation of the SS3 phase, which
undergoes a first order transition to the SF phase. Similar to the previous
cases discussed above, along with the gapless Goldstone mode, a gapped Higgs
excitation also exists in the SF phase, as depicted in Fig.12(a,b). When
approaching the SF-SS boundary, the lowest excitation branch in the SF phase
develops roton mode softening at $(k_{x},k_{y})=(\pi,\pi)$, $(\pi,0)$, and
$(0,\pi)$, with minimum roton gaps at $(\pi,0)$ and $(0,\pi)$, however, they
do not vanish at the phase boundary (see Fig.12(b)), reflecting the first
order nature of the SF-SS transition in this case. Unlike the case of SSCB-SF
continuous transition, the energy gap of the lowest gapped excitation in the
SS3 phase remains non-zero at the SS3-SF boundary as a consequence of such
first order transition. As shown in Fig.12(c), the low energy excitations of
SS3 phase consist of a gapless sound mode and a low energy gapped mode.
Finally, lowering the hopping strength leads to a stacked DW phase with a
configuration $(n_{0},0,0,0)$, exhibiting a gap opening at $\vec{k}=0$ in the
lowest excitation branch (see Fig.12(d)).

The above mentioned characteristics of the low-lying excitation spectrum in
the different regimes of interactions can be useful for the identification of
various phases and transition between them both at zero and finite
temperatures, which can also be experimentally relevant.

##  5 Summary of different supersolid phases and their stability

In this section, we summarize our results for the formation of various types
of supersolid (SS) phases of the hardcore as well as softcore bosons in a
square lattice. Below, we chart out the possible SS phases with different
density orderings and their regime of stability in Table1, as obtained within
mean field (MF) and cluster mean field (CMF) theory, which are also compared
with the available Quantum Monte Carlo (QMC) results. A schematic for the
different possible SS phases are shown in Fig.13.

Figure 13: Schematics of the different possible supersolid (SS) phases of
bosons in a square lattice. The circles with different sizes and shades of
gray color indicate non-integer fillings on the lattice sites.

Boson type | Interaction range | Methods | Supersolid (SS) phases
---|---|---|---
|  |  | ${\rm{\bf SS}}_{\rm{\bf CB}}$ | ${\rm{\bf SS}}_{\rm{\bf STR}}$ | SS3 | SS4
Hardcore bosons | $2V_{2}<V_{1}$ | MF | Stable | N.E | Stable | N.E
CMF | Stable***Although we find a region of stable checkerboard SS (${\rm SS}_{\rm CB}$) phase using $2\times 2$ CMF method, this region in phase diagram becomes smaller with increasing cluster size as more fluctuations are incorporated (see for example [54]). However, the QMC simulations reveal that, ${\rm SS}_{\rm CB}$ is unstable and phase separates into a mixture of checkerboard solid and superfluid phases [40]. | N.E | Stable | N.E
QMC | Unstable | N.E | Stable | N.E
$2V_{2}>V_{1}$ | MF | N.E | Stable | Stable | N.E
CMF | N.E | Stable | Stable | N.E
QMC | N.E | Stable | Stable | N.E
Softcore bosons | $2V_{2}<V_{1}$ | $U\lesssim 4V_{2}$ | MF | N.E | N.E | Stable | N.E
CMF | N.E | N.E | Stable | N.E
$4V_{2}\lesssim U\lesssim 4(V_{1}-V_{2})$ | MF | Stable | N.E | Stable | N.E
CMF | Stable | N.E | Stable | N.E
$U\gtrsim 4(V_{1}-V_{2})$ | MF | Stable | N.E | Stable | N.E
CMF | Stable | N.E | Stable | N.E
$2V_{2}>V_{1}$ | $U\lesssim 2V_{1}$ | MF | N.E | N.E | Stable | N.E
CMF | N.E | N.E | Stable | N.E
$2V_{1}\lesssim U\lesssim 4V_{2}$ | MF | N.E | Stable | Stable | Stable
CMF | N.E | Stable | Stable | Stable
$U\gtrsim 4V_{2}$ | MF | N.E | Stable | Stable | N.E†††SS4 does not exist at zero temperature in this regime, but it may exist as a stable phase at finite temperatures.
CMF | N.E | Stable | Stable | N.E† ‣ 1

Table 1: Summary of the stability of the various supersolid (SS) phases for
hardcore and softcore bosons in the different interaction regimes using mean
field (MF), cluster mean field (CMF) with $2\times 2$ cluster, and Quantum
Monte-Carlo (QMC) methods. Here, N.E stand for the phases which do not exist
as ground state in the specified interaction regime. Note that, the QMC
results are only available for certain cases.

Here we point out, even though the QMC simulations reveal that the
checkerboard SS phase (${\rm SS}_{\rm CB}$) remains unstable in the case of
hardcore bosons [40], it can be stabilized in the presence of a realistic long
range dipolar interaction [53]. In the present work, we have truncated the
long range interaction up to next nearest neighbor only, which is suitable for
numerical computations of CMF theory for clusters with finite size. However,
it has been shown in[54], within the CMF theory, that the checkerboard
ordering is favoured by considering an effective NN and NNN interaction
strength, which mimic the true long-range dipolar interaction. The long-range
nature of the dipolar interactions not only stabilize the checkerboard phases,
but also lead to the formation of various other density ordered phases at
lower fillings, exhibiting a devil’s staircase like structure [53]. In
addition, it can also give rise to different metastable configurations [47].
For larger lattice systems, the full long range dipolar interaction can also
be included using Ewald summation method, as done for QMC simulations in
Ref.[53, 56].

##  6 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 the 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
(CMF) theory to study the effect of quantum fluctuations systematically. We
consider bosons with nearest and next nearest neighbor interactions, which
play an essential role in the formation of the SS phase for softcore and
hardcore bosons respectively. Moreover, the interplay between nearest-neighbor
(NN), next nearest-neighbor (NNN), and on-site interactions lead to the
formation of different types of density ordering and structures in the SS
phase. We chart out various types of density ordering in solid and SS phases,
which can be categorized into three main cases. For sufficiently large on-site
repulsion, the competition between NN and NNN interactions gives rise to
checkerboard and striped density ordering. On the other hand, small on-site
repulsion favors a stacked structure with the accumulation of particles on a
single site of the unit cell. However, in an intermediate regime, various
types of solids can form with the coexistence of stacking and checkerboard (or
striped) order. For hardcore bosons, we study the finite temperature phases
only in the regime where NNN interaction dominates since the previous studies
reveal the instability of checkerboard supersolid, the analysis of which is
beyond the scope of the MF theory.

In the phase diagram, the SS phase appears in between the insulating density
wave (DW) and homogeneous superfluid (SF) phases. With increasing the
temperature, all the phases finally melt to the homogeneous normal fluid (NF)
phase. Melting of the SS phase is particularly interesting since it can occur
in at least two steps due to the presence of competing superfluidity and
density ordering. Different possible pathways for the transformation of the SS
to NF phase have been analyzed from the phase diagrams in the $T$-$t$ plane
for different regimes of interactions. In the $T$-$t$ plane of the phase
diagram, the boundaries of all the phases typically meet at a tetra-critical
point. However, from the cluster mean field results, we find that, for two
extreme limits of on-site interaction (for small $U$ as well in the hardcore
limit), 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 diagrams is similar to that obtained from the CMF
theory for intermediate range of $U$, the MF approximation overestimates the
critical temperatures, which can be reduced systematically by increasing the
cluster size. While it is expected that the quantum fluctuations will be
suppressed with increasing temperature, however, our results reveal that both
the thermal and quantum effects are important to understand the transitions
between various phases, as well as the stability of supersolids with different
orderings. In particular, the inclusion of correlations within CMF theory is
crucial to capture the nature of the tri/tetra-critical points at finite
temperatures. In addition, it would also be interesting to investigate the
modifications in the phase diagram that can arise due to the inclusion of the
long range nature of the dipolar interaction, which is however numerically
expensive to implement in the CMF theory, as it requires a larger cluster
size. Another possibility for studying such quantum phases at finite
temperature is to include the quadratic fluctuations on top of the simple mean
field solutions, which has recently been done for the Bose-Hubbard model using
a quantum Gutzwiller approach [98], and for hardcore bosons using the spin
wave theory [69].

We have also discussed how the different phases can be identified as well as
the transitions between them can be probed from the excitation spectrum, both
at zero and finite temperatures, which can be relevant for 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 the SS and DW phases, the excitation spectrum
of the SF phase develops a ‘roton’ minima, which signals possible breaking of
lattice translational symmetry. Moreover, the gap between the low-lying
excitation 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 the nature of density ordering.
At finite temperatures, we obtain the excitation spectrum from the
fluctuations of the equilibrium density matrix, which can also effectively
capture the melting process from the SS to NF phase. Apart from the gapless
sound mode, the appearance of a low-energy gapped excitation is an interesting
feature of the SS phase (analogous to the gapped mode in the dipolar
supersolids [34, 81]), and is reminiscent of the roton mode of the SF phase.
For the continuous transition from SSCB/SSSTR to the SF phase, both the gap of
the first excitation mode of the supersolid and the roton gap of the
superfluid vanish at the phase boundary.

In recent experiments, the signature of supersolidity has unambiguously been
observed in trapped dipolar quantum gas of ${}^{166}{\rm Er}$ and
${}^{164}{\rm Dy}$ atoms [19, 20, 21, 22, 23, 26, 27, 28, 16, 17, 18, 25, 24],
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 [20, 21, 25, 24]. The presence of a 2D optical lattice can stabilize
such supersolid phase of dipolar atoms. Moreover, tuning the on-site
interaction by Feshbach resonance can lead to different density ordering,
particularly the formation of a stacked SS phase. Also, the ultracold Rydberg
gas is another possible candidate for supersolids, due to their long range
interaction [15]. The 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
[99, 100]. Moreover, the tunability of effective interaction between Rydberg
dressed states is suitable for the formation of supersolids with different
density ordering. In recent experiments, the supersolid phase of dipolar gas
has been created by directly quenching the temperature, exhibiting sequential
formation of density order and superfluidity [21, 27, 28]. These experiments
show the importance to study the finite temperature phases of bosons and their
out-of-equilibrium dynamics. The pathways for the two step melting of the SS
phase 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 processes. The appearance of roton excitation as a
precursor of translational symmetry broken SS phase has been extensively
studied in experiments [16, 17, 18, 25, 30]. The translational symmetry broken
phases (such as DW and SS) and their melting at finite temperatures can also
be detected from the excitation branches, as discussed before.

In conclusion, we studied the finite temperature phases of hardcore as well as
softcore bosons using the mean field and cluster mean field theory, focusing
on the formation of supersolid phases with various density orderings. We have
demonstrated how the different phases as well as the transitions between them
at zero and finite temperatures can also be probed from the excitation
spectrum, which has importance for their detection in ongoing experiments.

## ACKNOWLEGMENTS

We thank S. Ray for stimulating and fruitful discussions.

##  Appendix A Finite cluster size scaling within Cluster mean field theory
and its accuracy

Here we check the quantitative accuracy of CMF theory at finite temperature by
performing the finite cluster-size scaling [54, 77]. We consider the scaling
parameter, $\lambda=\frac{N_{B}}{N_{\mathcal{C}}\times z/2}$ (introduced in
Ref.[77]), with $N_{B}$ and $N_{\mathcal{C}}$ denoting the number of bonds and
lattice sites within the cluster $\mathcal{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\times 2$, $2\times 3$, $2\times 4$ (see
Fig.14(a) for schematics), and finally extrapolate the values at the
thermodynamic limit $\lambda\to 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 CMF theory at finite temperature using finite cluster-size
scaling, with the Quantum Monte-Carlo (QMC) results of Ref.[101].

Figure 14: (a) Schematics of the different cluster sizes in 2D. (b,c) Finite
cluster-size scaling analysis to capture the critical temperatures of DW-NF
transition for different set of parameters with $V_{2}=0$, by using different
cluster sizes. The solid red line is the linear fit which gives the
extrapolated value at $\lambda\to 1$ under CMF theory. The Quantum Monte-Carlo
(QMC) values are marked by a red cross in (b,c), and are taken from Ref.[101].
The extrapolated CMF (QMC) values in (b) and (c) are $T^{\rm DW-
NF}_{c}/V_{1}\simeq 0.664$ $(0.567)$ and $T^{\rm DW-NF}_{c}/t\simeq 1.785$
$(1.5)$, respectively.

As shown in Fig.14(b,c), the transition points acquired from CMF theory, after
performing the linear extrapolation up to $\lambda\to 1$, differ significantly
from the single site MF results and approach the more accurate values obtained
from QMC. As evident from above, with the larger cluster sizes like $4\times
4$ and so on, the CMF theory can yield results having a better agreement with
the QMC values, which can improve the phase boundaries obtained from MF
theory. 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 [97], however the transition point obtained from
extrapolation after using the finite cluster size scaling approaches the QMC
value.

Figure 15: Variation of (a) density imbalance $\delta_{12}$ capturing the DW-
NF transition and (b) SF order $\bar{\phi}$ capturing the SF-NF transition in
3D using a $2\times 2\times 2$ cubic cluster (filled circles) and mean field
approximation (dotted line). Parameters chosen: $V_{1}/U=1/6$, $\mu/U=0.7$,
and $V_{2}/U=0$. The QMC values marked by a red cross in (a,b) are taken from
Ref.[49]. The CMF (QMC) values in (a) and (b) are $T^{\rm DW-NF}_{c}/U\simeq
0.26$ $(0.24)$ and $T^{\rm SF-NF}_{c}/U\simeq 0.383$ $(0.35)$.

The CMF theory 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\times 2\times 2$ cubic cluster, and compare the critical temperatures with
that of MF and QMC results of Ref.[49]. To implement CMF method for bosons
with longer range interaction, larger cluster size can be required for
extracting the transition points more accurately. In Fig.15(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 CMF theory improves
the MF value and are much closer to the QMC results. In three dimensions, due
to the reduced effect of fluctuations, the finite temperature CMF theory
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.

##  Appendix B Excitation spectrum of hardcore bosons

In this appendix, we provide the (semi) analytical expressions of low-lying
collective excitations for different phases of hardcore bosons in the presence
of NN and NNN interactions at zero (finite) temperatures.

###  B.1 Zero temperature

Following the method outlined in Sec.2, here we provide the detailed
calculations and derive the analytical expressions of excitation spectrum
$\omega(\vec{k})$ for different phases of hardcore bosons.

The excitation branches at zero temperature can be calculated analytically
from the equations of motion (EOM) of the time dependent Gutzwiller
amplitudes, given by Eq.(5). For hardcore bosons with $n_{i}\in\\{0,1\\}$, the
EOM for $f^{(1)}_{i}$ and $f^{(0)}_{i}$ at $i^{th}$ site are given by,

$\displaystyle\imath\dot{f}^{(0)}_{i}=-tf^{(1)}_{i}(\sum^{\rm
NN}_{j}f^{*(1)}_{j}f^{(0)}_{j})-\lambda_{i}f^{(0)}_{i}$ (16a)
$\displaystyle\imath\dot{f}^{(1)}_{i}=-tf_{i}^{(0)}(\sum^{\rm
NN}_{j}f^{*(0)}_{j}f^{(1)}_{j})-(\mu+\lambda_{i})f_{i}^{(1)}\\!+\\!V_{1}f_{i}^{(1)}\\!\\!\sum^{\rm
NN}_{j}|f_{j}^{(1)}|^{2}\\!+\\!V_{2}f_{i}^{(1)}\sum^{\rm
NNN}_{j^{\prime}}|f_{j^{\prime}}^{(1)}|^{2}$ (16b)

Below, we provide the analytical expressions for dispersion of different
phases of hardcore bosons at zero temperature and finite hopping strength
$t\neq 0$.

Checkerboard DW-1/2: By considering the two sub-lattice density structure
$(n_{A},n_{B})\equiv(1,0)$ with checkerboard ordering (cf. Fig.1(a)), the
steady state values of Gutzwiller amplitudes are $\bar{f}_{A}^{(0)}=0$,
$\bar{f}_{A}^{(1)}=1$ and $\bar{f}_{B}^{(0)}=1$, $\bar{f}_{B}^{(1)}=0$. By
substituting the steady state values in Eq.(16) and using the steady state
condition, we obtain the Lagrange multipliers at each sub-lattice, given by
$\lambda_{A}=-\mu+4V_{2}$ and $\lambda_{B}=0$. After introducing the small
amplitude oscillations around the steady state values $\bar{f}^{(n)}_{i}$(see
Eq.(6)), the linearized fluctuation equations in the momentum space within the
reduced Brillouin zone can be written as following,

$\displaystyle[\,\omega-(\mu-4V_{2})\,]\delta
f_{A}^{(0)}\\!=\\!-\epsilon(\vec{k})\delta f_{B}^{*(1)}$ (17a)
$\displaystyle[\,\omega+(\mu-4V_{1})\,]\delta
f_{B}^{(1)}\\!=\\!-\epsilon(\vec{k})\delta f_{A}^{*(0)}$ (17b)
$\displaystyle[\,\omega+(\mu-4V_{2})\,]\delta
f_{A}^{*(0)}\\!=\\!\epsilon(\vec{k})\delta f_{B}^{(1)}$ (17c)
$\displaystyle[\,\omega-(\mu-4V_{1})\,]\delta
f_{B}^{*(1)}\\!=\\!\epsilon(\vec{k})\delta f_{A}^{(0)}$ (17d)

where $\epsilon(\vec{k})=2t(\cos k_{x}+\cos k_{y})$. Solving the above
equations simultaneously for $\omega(\vec{k})$ yields the particle ($+$) and
hole ($-$) excitation $\omega_{\pm}(\vec{k})$ for checkerboard DW phase with
filling $n_{0}=1/2$, which is given by,

$\displaystyle\omega_{\pm}(\vec{k})=\pm\tilde{\mu}\\!+\\!\sqrt{4(V_{1}-V_{2})^{2}-\epsilon(\vec{k})^{2}}$
(18)

where $\tilde{\mu}=-\mu+2(V_{1}+V_{2})$.

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.(16), the Lagrange multipliers at $A$
and $B$ sub-lattice can be fixed as $\lambda_{A}=-\mu+2V_{1}$ and
$\lambda_{B}=0$, respectively. The linearized fluctuation equations in
momentum space within the reduced Brillouin zone can be written as following,

$\displaystyle[\,\omega+(\epsilon_{2}(\vec{k})\\!+\\!2V_{1}\\!-\\!\mu)\,]\delta
f_{A}^{(0)}\\!=\\!-\epsilon_{1}(\vec{k})\delta f_{B}^{*(1)}$ (19a)
$\displaystyle[\,\omega+(\epsilon_{2}(\vec{k})\\!-\\!4V_{2}\\!-\\!2V_{1}\\!+\\!\mu)\,]\delta
f_{B}^{(1)}\\!=\\!-\epsilon_{1}(\vec{k})\delta f_{A}^{*(0)}$ (19b)
$\displaystyle[\,\omega-(\epsilon_{2}(\vec{k})\\!+\\!2V_{1}\\!-\\!\mu)\,]\delta
f_{A}^{*(0)}\\!=\\!\epsilon_{1}(\vec{k})\delta f_{B}^{(1)}$ (19c)
$\displaystyle[\,\omega-(\epsilon_{2}(\vec{k})\\!-\\!4V_{2}\\!-\\!2V_{1}\\!+\\!\mu)\,]\delta
f_{B}^{*(1)}\\!=\\!\epsilon_{1}(\vec{k})\delta f_{A}^{(0)}$ (19d)

where $\epsilon_{1}(\vec{k})=2t\cos k_{x}$, and $\epsilon_{2}(\vec{k})=2t\cos
k_{y}$. By solving the above equations simultaneously for $\omega(k)$, the
particle ($+$) and hole ($-$) excitation $\omega_{\pm}(\vec{k})$ for striped
DW phase with filling $n_{0}=1/2$ is given by,

$\displaystyle\omega_{\pm}(\vec{k})=\pm\tilde{\mu}\\!+\\!\sqrt{\epsilon_{2}(\vec{k})^{2}-\epsilon_{1}(\vec{k})^{2}-4V_{2}(\epsilon_{2}(\vec{k})-V_{2})}$
(20)

Homogeneous SF phase: For a homogeneous SF phase, the Gutzwiller variational
wavefunction at each site can be written as
$|\psi\rangle\\!=\\!f^{(0)}|0\rangle\\!+\\!f^{(1)}|1\rangle$, where we can use
the parametrization $f^{(0)}\\!=\\!\sin{(\theta/2)}$,
$f^{(1)}\\!=\\!\cos{(\theta/2)}$. From the steady state conditions of the EOM,
we obtain $\lambda=-4t\cos^{2}{(\theta/2)}$ with
$\theta=\cos^{-1}{\left(\frac{\mu-2(V_{1}+V_{2})}{2(V_{1}+V_{2}+2t)}\right)}$,
and subsequently solve the following linearized fluctuation equations in the
momentum space within the extended Brillouin zone,

$\displaystyle\left[\omega\\!+\\!\tilde{\epsilon}(\vec{k})\cos^{2}{\\!\frac{\theta}{2}}\right]\\!\frac{\delta\\!f^{(0)}}{\sin{\theta}}\\!=\\!-\\!2t\delta
f^{(1)}\\!-\\!\frac{\epsilon(\vec{k})}{2}\delta f^{*(1)}$ (21)
$\displaystyle\left[\omega\\!-\\!\tilde{\epsilon}(\vec{k})\cos^{2}{\\!\frac{\theta}{2}}\right]\\!\frac{\delta\\!f^{*(0)}}{\sin{\theta}}=\\!2t\delta
f^{*(1)}\\!+\\!\frac{\epsilon(\vec{k})}{2}\delta f^{(1)}$ (22)
$\displaystyle\left[\omega\\!+\\!\mathcal{E}(\vec{k})\right]\\!\frac{\delta
f^{(1)}}{\sin{\theta}}\\!-\\!X(\vec{k})\\!\cos^{2}{\\!\frac{\theta}{2}}\frac{\delta
f^{*(1)}}{\sin{\theta}}\\!=-\\!2t\delta
f^{(0)}\\!-\\!\frac{\epsilon(\vec{k})}{2}\delta f^{*(0)}$ (23)
$\displaystyle\left[\omega\\!-\\!\mathcal{E}(\vec{k})\right]\\!\frac{\delta
f^{*(1)}}{\sin{\theta}}\\!+\\!X(\vec{k})\\!\cos^{2}{\\!\frac{\theta}{2}}\frac{\delta
f^{(1)}}{\sin{\theta}}\\!=\\!2t\delta
f^{*(0)}\\!+\\!\frac{\epsilon(\vec{k})}{2}\delta f^{(0)}$ (24)

where $\tilde{\epsilon}(\vec{k})=\epsilon(\vec{k})-4t$,
$\mathcal{E}({\vec{k}})=\epsilon(\vec{k})+\mu-[\epsilon(\vec{k})+4t+X(\vec{k})+4(V_{1}+V_{2})]\cos^{2}{(\theta/2)}$,
and $X(\vec{k})=V_{1}(\cos k_{x}+\cos
k_{y})+V_{2}\left[\cos(k_{x}+k_{y})+\cos(k_{x}-k_{y})\right]$. The excitation
spectrum $\omega(\vec{k})$ is given by,

$\displaystyle\omega(\vec{k})=\sqrt{\cos^{2}{\theta}\,\tilde{\epsilon}(\vec{k})^{2}-\sin^{2}{\theta}\,\tilde{\epsilon}(\vec{k})\\!\left(4t+X(\vec{k})\right)}$
(25)

###  B.2 Finite temperature

The method of linear stability analysis can be extended to finite
temperatures, where the stability analysis of the finite temperature density
matrix $\hat{\rho}_{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,

$\displaystyle\hat{\rho}_{i}=\begin{bmatrix}\left(1-m_{i}\right)/2&\phi^{*}_{i}\\\
\phi_{i}&\left(1+m_{i}\right)/2\end{bmatrix}$ (26)

using the above density matrix, the SF order parameter is given by,
$\phi_{i}={\rm Tr}(\hat{\rho}_{i}\hat{a}_{i})$ and the average density can be
written in terms of $m_{i}$ as, $\langle\hat{n}_{i}\rangle={\rm
Tr}(\hat{\rho}_{i}\hat{n_{i}})=(1+m_{i})/2$. By introducing a small amplitude
fluctuation $\delta\hat{\rho}_{i}(t)$ around the equilibrium density matrix
$\hat{\rho}^{eq}_{i}$, such that
$\hat{\rho}_{i}(t)=\hat{\rho}^{eq}_{i}+\delta\hat{\rho}_{i}(t)$, and plugging
it into Eq.(8), we obtain the following,

$\displaystyle\delta\dot{\hat{\rho}}_{i}=-\imath[\hat{\mathcal{H}}^{\rm
MF}_{i},\delta\hat{\rho}_{i}]-\imath[\delta\hat{\mathcal{H}}^{\rm
MF}_{i},\hat{\rho}^{eq}_{i}]$ (27)

where $\delta\hat{\mathcal{H}}^{\rm MF}_{i}$ contains the fluctuations around
the mean field terms in $\hat{\mathcal{H}}^{\rm MF}_{i}$. Substituting
$\delta\hat{\rho}_{i}(t)=e^{\imath(\vec{k}.\vec{r}_{i}-\omega(\vec{k})t)}\delta\hat{\rho}(\vec{k})$
in Eq.(27) and retaining terms only up to the first order of
$\delta\hat{\rho}(\vec{k})$, we obtain the fluctuation matrix corresponding to
the set of linearized equations in momentum space. The real part of the
eigenvalues $\omega(\vec{k})$ of the fluctuation matrix provides the
excitation branches, whereas ${\rm Im}[\omega(\vec{k})]=0$ ensures the
stability of $\hat{\rho}^{eq}$. Below, we provide the semi analytical
expressions for low-lying excitations of some of the phases of hardcore bosons
at finite temperature.

Checkerboard DW: The checkerboard DW (DWCB) is an inhomogeneous phase with 2
sub-lattice structure and checkerboard ordering (cf. Fig.1(a)). In order to
obtain the excitation spectrum $\omega(\vec{k})$, it requires solving a set of
six fluctuation equations. The low-lying excitation branches $\omega_{\pm}(k)$
for DWCB at finite temperature is given by the following,

$\displaystyle\omega_{\pm}(\vec{k})\\!=\\!\frac{\pm(\Gamma_{A}\\!+\\!\Gamma_{B})\\!+\\!\sqrt{(\Gamma_{A}\\!-\\!\Gamma_{B})^{2}\\!+\\!4\mathcal{K}_{A}\mathcal{K}_{B}}}{2}$
(28)

where $\Gamma_{i}(\vec{k})=2V_{1}(1+m_{\bar{i}})+2V_{2}(1+m_{i})-\mu$,
$\mathcal{K}_{i}(\vec{k})=\epsilon(\vec{k})m_{i}$, and
$\epsilon(\vec{k})=2t(\cos k_{x}+\cos k_{y})$, for $i,\,\bar{i}=A,B$ with
$i\neq\bar{i}$. Note that, here the parameter $m_{i}$ has an implicit
dependence on temperature, which can be obtained numerically in a self
consistent manner. It is easy to check, by setting $m_{A}=1$ and $m_{B}=-1$
for the DWCB phase in Eq.28, the expression matches with that in Eq.(18) for
two sub-lattice density structure $(n_{A},n_{B})=(1,0)$ with filling
$n_{0}=1/2$ at zero temperature. Note that, one can also obtain the excitation
spectrum for checkerboard SS (SSCB) at finite temperature within MF
approximation, however as mentioned in the main text, it has been known from
previous QMC studies that, the SSCB phase is unstable even in the presence of
NNN interactions [40].

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 linearized fluctuation
equations to obtain the low-lying excitations $\omega_{\pm}(\vec{k})$ for
stripe SS phase at finite temperature, which is given by the following,

$\displaystyle\omega_{\pm}(\vec{k})=\left[\frac{\mathcal{A}^{\prime}(\vec{k})\pm\sqrt{\mathcal{B}^{\prime}(\vec{k})+\mathcal{C}^{\prime}(\vec{k})+\mathcal{D}^{\prime}(\vec{k})}}{2}\right]^{1/2}$
(29)

where $\mathcal{A}^{\prime}$, $\mathcal{B}^{\prime}$, $\mathcal{C}^{\prime}$,
$\mathcal{D}^{\prime}$ are given by,

$\displaystyle\mathcal{A}^{\prime}=\,\,$
$\displaystyle\Gamma^{{}^{\prime}2}_{A}\\!+\\!\Gamma^{{}^{\prime}2}_{B}\\!+\\!2(\beta^{{}^{\prime}}_{B}\zeta^{{}^{\prime}}_{A}\\!+\\!\alpha^{{}^{\prime}}_{A}\eta^{{}^{\prime}}_{A}\\!+\\!\alpha^{{}^{\prime}}_{B}\eta^{{}^{\prime}}_{B}\\!+\\!\mathcal{K}^{{}^{\prime}}_{A}\mathcal{K}^{{}^{\prime}}_{B}\\!+\\!\beta^{{}^{\prime}}_{A}\chi^{{}^{\prime}}_{B})$
(30a) $\displaystyle\mathcal{B}^{\prime}=\,\,$
$\displaystyle\Gamma^{{}^{\prime}4}_{A}\\!+\\!\Gamma^{{}^{\prime}4}_{B}\\!-\\!4\Gamma^{{}^{\prime}2}_{B}(\beta^{{}^{\prime}}_{B}\zeta^{{}^{\prime}}_{A}+\alpha^{{}^{\prime}}_{A}\eta^{{}^{\prime}}_{A}-\alpha^{{}^{\prime}}_{B}\eta^{{}^{\prime}}_{B}-\mathcal{K}^{{}^{\prime}}_{A}\mathcal{K}^{{}^{\prime}}_{B}$
$\displaystyle\\!-\\!\beta^{{}^{\prime}}_{A}\chi^{{}^{\prime}}_{B})+8[\Gamma^{{}^{\prime}}_{A}\Gamma^{{}^{\prime}}_{B}\mathcal{K}^{{}^{\prime}}_{A}\mathcal{K}^{{}^{\prime}}_{B}+(\Gamma^{{}^{\prime}}_{B}\\!+\\!\Gamma^{{}^{\prime}}_{A})(\beta^{{}^{\prime}}_{B}\eta^{{}^{\prime}}_{B}\mathcal{K}^{{}^{\prime}}_{A}$
$\displaystyle\\!+\\!\alpha^{{}^{\prime}}_{B}\zeta^{{}^{\prime}}_{A}\mathcal{K}^{{}^{\prime}}_{B}\\!+\\!\beta^{{}^{\prime}}_{A}\eta^{{}^{\prime}}_{A}\mathcal{K}^{{}^{\prime}}_{B}\\!+\\!\alpha^{{}^{\prime}}_{A}\mathcal{K}^{{}^{\prime}}_{A}\chi^{{}^{\prime}}_{B})]$
(30b) $\displaystyle\mathcal{C}^{\prime}=\,\,$ $\displaystyle
4(\beta^{{}^{\prime}2}_{B}\zeta^{{}^{\prime}2}_{A}+\beta^{{}^{\prime}2}_{A}\chi^{{}^{\prime}2}_{B})-2\Gamma^{{}^{\prime}2}_{A}[\Gamma^{{}^{\prime}2}_{B}\\!-\\!2(\beta^{{}^{\prime}}_{B}\zeta^{{}^{\prime}}_{A}\\!+\\!\alpha^{{}^{\prime}}_{A}\eta^{{}^{\prime}}_{A}$
$\displaystyle\\!-\\!\alpha^{{}^{\prime}}_{B}\eta^{{}^{\prime}}_{B}\\!+\\!\mathcal{K}^{{}^{\prime}}_{A}\mathcal{K}^{{}^{\prime}}_{B}-\beta^{{}^{\prime}}_{A}\chi^{{}^{\prime}}_{B})]\\!+\\!4(\alpha^{{}^{\prime}}_{A}\eta^{{}^{\prime}}_{A}\\!-\\!\alpha^{{}^{\prime}}_{B}\eta^{{}^{\prime}}_{B})^{2}$
(30c) $\displaystyle\mathcal{D}^{\prime}=\,\,$ $\displaystyle
8(2\alpha^{{}^{\prime}}_{A}\alpha^{{}^{\prime}}_{B}\zeta^{{}^{\prime}}_{A}\\!+\\!\alpha^{{}^{\prime}}_{A}\beta^{{}^{\prime}}_{A}\eta^{{}^{\prime}}_{A}\\!+\\!\alpha^{{}^{\prime}}_{B}\beta^{{}^{\prime}}_{A}\eta^{{}^{\prime}}_{B})\chi^{{}^{\prime}}_{B}\\!+\\!8\beta^{{}^{\prime}}_{B}$
$\displaystyle(2\beta^{{}^{\prime}}_{A}\eta^{{}^{\prime}}_{A}\eta^{{}^{\prime}}_{B}\\!+\\!\alpha^{{}^{\prime}}_{A}\zeta^{{}^{\prime}}_{A}\eta^{{}^{\prime}}_{A}\\!+\\!\alpha^{{}^{\prime}}_{B}\zeta^{{}^{\prime}}_{A}\eta^{{}^{\prime}}_{B}-\beta^{{}^{\prime}}_{A}\zeta^{{}^{\prime}}_{A}\chi^{{}^{\prime}}_{B})$
(30d)

with
$\alpha^{{}^{\prime}}_{i}(\vec{k})=2t[2(\phi_{i}+\phi_{\bar{i}})-\phi_{i}\epsilon_{2}(\vec{k})/t]$,
$\beta^{{}^{\prime}}_{i}(\vec{k})=-2\epsilon_{1}(\vec{k})\phi_{i}$,
$\zeta^{{}^{\prime}}_{i}(\vec{k})=\phi_{i}(X(\vec{k})-V_{1}\cos{k_{y}})$,
$\chi^{{}^{\prime}}_{i}(\vec{k})=\phi_{i}V_{2}[\cos{(k_{x}+k_{y})}+\cos{(k_{x}-k_{y})}]+\phi_{i}V_{1}\epsilon_{1}(\vec{k})/2t$,
$\mathcal{K}^{{}^{\prime}}_{i}(\vec{k})=\epsilon_{1}(\vec{k})m_{i}$,
$\Gamma^{{}^{\prime}}_{i}(\vec{k})=\epsilon_{2}(\vec{k})m_{i}+V_{1}(2+m_{i}+m_{\bar{i}})+2V_{2}(1+m_{\bar{i}})-\mu$,
$\eta^{{}^{\prime}}_{i}(\vec{k})=2t(\phi_{i}+\phi_{\bar{i}})+V_{1}\phi_{i}\epsilon_{2}(\vec{k})/2t$,
$\epsilon_{1}(\vec{k})=2t\cos k_{x}$, and $\epsilon_{2}(\vec{k})=2t\cos
k_{y}$, for $i=A,B$ with $i\neq\bar{i}$. Here, $\phi_{i}$ and $m_{i}$ at
finite temperature can be obtained self consistently from numerics. To get the
excitation spectrum for stripe DW (DWSTR) phase, we can substitute
$\phi_{i}=0$ in the above equations, which gives the following,

$\displaystyle\omega_{\pm}(\vec{k})\\!=\\!\frac{\pm({\Gamma}^{{}^{\prime}}_{A}\\!+\Gamma^{{}^{\prime}}_{B})\\!+\\!\sqrt{(\Gamma^{{}^{\prime}}_{A}\\!-\\!\Gamma^{{}^{\prime}}_{B})^{2}\\!+\\!4\mathcal{K}^{{}^{\prime}}_{A}\mathcal{K}^{{}^{\prime}}_{B}}}{2}$
(31)

Note that, by setting $m_{A}=1$ and $m_{B}=-1$ for DWSTR phase, the above
expression matches with that in Eq.(20), for two sub-lattice density structure
$(n_{A},n_{B})=(1,0)$ with filling $n_{0}=1/2$ at zero temperature.

Homogeneous SF phase: Within the homogeneous SF phase, the set of linear
fluctuation equations in the momentum space can be written as following,

$\displaystyle\omega\,\delta m=2\tau(\vec{k})[\,\delta\phi-\delta\phi^{*}\,]$
(32a)
$\displaystyle[\,\omega-2\gamma(\vec{k})\,]\delta\phi=2\Lambda(\vec{k})\delta
m$ (32b)
$\displaystyle[\,\omega+2\gamma(\vec{k})\,]\delta\phi^{*}=-2\Lambda(\vec{k})\delta
m$ (32c)

where $\tau(\vec{k})=-\phi\tilde{\epsilon}(\vec{k})$,
$\Lambda(\vec{k})=\phi(2t+X(\vec{k})/2)$,
$\gamma(\vec{k})=m\epsilon(\vec{k})/2-\mu/2+(V_{1}+V_{2})(1+m)$, and
$\tilde{\epsilon}(\vec{k})=\epsilon(\vec{k})-4t$. Similar to the previous
phases, the parameters $\phi$ and $m$ can be obtained self consistently. By
solving the set of fluctuation equations in Eq.(32), we obtain the lowest
lying excitation branch for homogeneous SF phase at finite temperature,

$\displaystyle\omega({\vec{k}})=2\sqrt{2\tau(\vec{k})\Lambda(\vec{k})+\gamma(\vec{k})^{2}}.$
(33)

One can easily check, at zero temperature, when
$m=\cos{\theta}=\left(\frac{\mu-2(V_{1}+V_{2})}{2(V_{1}+V_{2}+2t)}\right)$ and
$\phi=(\sin{\theta})/2$, the above expression in Eq.(33) agrees with that in
Eq.(25).

## References

## References

  * [1] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature (London) 415, 39 (2002).
  * [2] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen, and U. Sen, Adv. Phys. 56, 243 (2007).
  * [3] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * [4] T. Langen, R. Geiger, and J. Schmiedmayer, Annu. Rev. Condens. Matter Phys. 6, 201 (2015).
  * [5] O. Penrose and L. Onsager, Phys. Rev. 104, 576 (1956).
  * [6] A. F. Andreev and I. M. Lifshitz, Sov. Phys. JETP 29, 1107 (1969).
  * [7] G. V. Chester, Phys. Rev. A 2, 256 (1970).
  * [8] A. J. Leggett, Phys. Rev. Lett. 25, 1543 (1970).
  * [9] K.-S. Liu and M. E. Fisher, J. Low Temp. Phys. 10, 655 (1972).
  * [10] M. Boninsegni and N. V. Prokof’ev, Rev. Mod. Phys. 84, 759 (2012).
  * [11] E. Kim and M. H. W. Chan, Nature (London) 427, 225 (2004); D. Y. Kim and M. H. W. Chan, Phys. Rev. Lett. 109, 155301 (2012).
  * [12] J. Stuhler, A. Griesmaier, T. Koch, M. Fattori, T. Pfau, S. Giovanazzi, P. Pedri, and L. Santos, Phys. Rev. Lett. 95, 150406 (2005).
  * [13] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, Rep. Prog. Phys. 72, 126401 (2009).
  * [14] L. Chomaz, I. Ferrier-Barbut, F. Ferlaino, B. Laburthe-Tolra, B. Lev, and T. Pfau, Rep. Prog. Phys. 86, 026401 (2022).
  * [15] R. Löw, H. Weimer, J. Nipper, J. B. Balewski, B. Butscher, H. P. Büchler, and T. Pfau, J. Phys. B 45, 113001 (2012).
  * [16] L. Chomaz, R. M. W. van Bijnen, D. Petter, G. Faraoni, S. Baier, J. H. Becher, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, Nat. Phys. 14, 442 (2018).
  * [17] D. Petter, G. Natale, R. M. W. van Bijnen, A. Patscheider, M. J. Mark, L. Chomaz, and F. Ferlaino, Phys. Rev. Lett. 122, 183401 (2019).
  * [18] G. Natale, R. M. W. van Bijnen, A. Patscheider, D. Petter, M. J. Mark, L. Chomaz, and F. Ferlaino, Phys. Rev. Lett. 123, 050402 (2019).
  * [19] L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, Phys. Rev. Lett. 122, 130405 (2019).
  * [20] F. Böttcher, J. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, Phys. Rev. X 9, 011051 (2019).
  * [21] L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen, M. J. Mark, F. Ferlaino, Phys. Rev. X 9, 021012 (2019).
  * [22] L. Tanzi, S. M. Roccuzzo, E. Lucioni, F. Famà, A. Fioretti, C. Gabbanini, G. Modugno, A. Recati, and S. Stringari, Nature (London) 574, 382 (2019).
  * [23] M. Guo, F. Böttcher, J. Hertkorn, J.-N. Schmidt, M. Wenzel, H. P. Büchler, T. Langen, and T. Pfau, Nature (London) 574, 386 (2019).
  * [24] P. Ilzhöfer, M. Sohmen, G. Durastante, C. Politi, A. Trautmann, G. Natale, G. Morpurgo, T. Giamarchi, L. Chomaz, M. J. Mark, and F. Ferlaino, Nat. Phys. 17, 356 (2021).
  * [25] J. Hertkorn, J.-N. Schmidt, F. Böttcher, M. Guo, M. Schmidt, K. S. H. Ng, S. D. Graham, H. P. Büchler, T. Langen, M. Zwierlein, and T. Pfau, Phys. Rev. X 11, 011037 (2021).
  * [26] M. A. Norcia, C. Politi, L. Klaus, E. Poli, M. Sohmen, M. J. Mark, R. N. Bisset, L. Santos, and F. Ferlaino, Nature (London) 596, 357 (2021).
  * [27] M. Sohmen, C. Politi, L. Klaus, L. Chomaz, M. J. Mark, M. A. Norcia, and F. Ferlaino, Phys. Rev. Lett. 126, 233401 (2021).
  * [28] T. Bland, E. Poli, C. Politi, L. Klaus, M. A. Norcia, F. Ferlaino, L. Santos, R. N. Bisset, Phys. Rev. Lett. 128, 195302 (2022).
  * [29] K. Baumann, C. Guerlin, F. Brennecke, and T. Esslinger, Nature (London) 464, 1301 (2010).
  * [30] R. Mottl, F. Brennecke, K. Baumann, R. Landig, T. Donner, and T. Esslinger, Science 336, 1570 (2012).
  * [31] J. Klinder, H. Keßler, M. Wolke, L. Mathey, and A. Hemmerich, Proc. Natl. Acad. Sci. U.S.A. 112, 3290 (2015).
  * [32] R. Landig, L. Hruby, N. Dogra, M. Landini, R. Mottl, T. Donner, and T. Esslinger, Nature (London) 532, 476 (2016).
  * [33] J. Léonard, A. Morales, P. Zupancic, T. Esslinger, and T. Donner, Nature (London) 543, 87 (2017).
  * [34] J. Léonard, A. Morales, P. Zupancic, T. Donner, and T. Esslinger, Science 358, 1415 (2017).
  * [35] J. R. Li, J. Lee, W. Huang, S. Burchesky, B. Shteynas, F. Ç. Topi, A. O. Jamison, and W. Ketterle, Nature (London) 543, 91 (2017).
  * [36] Y.-J. Lin, K. Jiménez-García, and I. B. Spielman, Nature (London) 471, 83 (2011).
  * [37] A. van Otterlo and K. H. Wagenblast, Phys. Rev. Lett. 72, 3598 (1994).
  * [38] A. van Otterlo, K. H. Wagenblast, R. Baltin, C. Bruder, R. Fazio, and G. Schön, Phys. Rev. B 52, 16176 (1995).
  * [39] G. G. Batrouni, R. T. Scalettar, G. T. Zimanyi, and A. P. Kampf, Phys. Rev. Lett. 74, 2527 (1995); R. T. Scalettar, G. G. Batrouni, A. P. Kampf, and G. T. Zimanyi, Phys. Rev. B 51, 8467 (1995).
  * [40] G. G. Batrouni and R. T. Scalettar, Phys. Rev. Lett. 84, 1599 (2000); F. Hébert, G. G. Batrouni, R. T. Scalettar, G. Schmid, M. Troyer, and A. Dorneich, Phys. Rev. B 65, 014513 (2001).
  * [41] C. Pich and E. Frey, Phys. Rev. B 57, 13712 (1998).
  * [42] K. Góral, L. Santos, and M. Lewenstein, Phys. Rev. Lett. 88, 17 (2002); S. Yi, T. Li, and C. P. Sun, Phys. Rev. Lett. 98, 260405 (2007).
  * [43] G. Schmid and M. Troyer, Phys. Rev. Lett. 93, 067003 (2004).
  * [44] P. Sengupta, L. P. Pryadko, F. Alet, M. Troyer, and G. Schmid, Phys. Rev. Lett. 94, 207202 (2005).
  * [45] D. L. Kovrizhin, G. V. Pai, and S. Sinha, Europhys. Lett. 72, 162 (2005).
  * [46] N. Laflorencie and F. Mila, Phys. Rev. Lett. 99, 027202 (2007).
  * [47] C. Menotti, C. Trefzger, and M. Lewenstein, Phys. Rev. Lett. 98, 235301 (2007).
  * [48] Y. Chen, R. G. Melko, S. Wessel, and Y. Kao, Phys. Rev. B 77, 014524 (2008).
  * [49] K. Yamamoto, S. Todo, and S. Miyashita, Phys. Rev. B 79, 094503 (2009).
  * [50] K.-K. Ng, Phys. Rev. B 82, 184505 (2010).
  * [51] A. Hubener, M. Snoek, and W. Hofstetter, Phys. Rev. B 80, 245109 (2009).
  * [52] L. Mathey, I. Danshita, C. W. Clark, Phys. Rev. A 79, 011602(R) (2009).
  * [53] B. Capogrosso-Sansone, C. Trefzger, M. Lewenstein, P. Zoller, and G. Pupillo, Phys. Rev. Lett. 104, 125301 (2010).
  * [54] D. Yamamoto, A. Masaki, and I. Danshita, Phys. Rev. B 86, 054516 (2012).
  * [55] T. Ohgoe, T. Suzuki, and N. Kawashima, Phys. Rev. B 86, 054520 (2012).
  * [56] D. Grimmer, A. Safavi-Naini, B. Capogrosso-Sansone, and Ş. G. Söyler, Phys. Rev. A 90, 043635 (2014).
  * [57] K. Saha, S. Sinha, and K. Sengupta, Phys. Rev. A 89, 023618 (2014).
  * [58] F. Heydarinasab and J. Abouie, Phys. Rev. B 96, 104406 (2017).
  * [59] O. Nguyen and L. Dang, Eur. Phys. J. B 90, 71 (2017).
  * [60] S. Bandyopadhyay, R. Bai, S. Pal, K. Suthar, R. Nath, and D. Angom, Phys. Rev. A 100, 053623 (2019).
  * [61] K. Suthar, H. Sable, R. Bai, S. Bandyopadhyay, S. Pal, and D. Angom, Phys. Rev. A 102, 013320 (2019).
  * [62] K. Suthar, R. Kraus, H. Sable, D. Angom, G. Morigi, and J. Zakrzewski, Phys. Rev. B 102, 214503 (2020).
  * [63] G. Masella, A. Angelone, F. Mezzacapo, G. Pupillo, and N. V. Prokof’ev, Phys. Rev. Lett. 123, 045301 (2019); A. Angelone, T. Ying, F. Mezzacapo, G. Masella, M. Dalmonte, and G. Pupillo, Phys. Rev. A 101, 063603 (2020).
  * [64] S. Prestipino, Phys. Rev. A 103, 033313 (2021).
  * [65] Y. Li, A. Geißler, W. Hofstetter, and W. Li, Phys. Rev. A 97, 023619 (2018); M. Barbier, H. Lütjeharms, and W. Hofstetter, Phys. Rev. A 105, 013326 (2022).
  * [66] T. D. Kühner and H. Monien, Phys. Rev. B 58, R14741(R) (1998); T. D. Kühner, S. R. White, and H. Monien, Phys. Rev. B 61, 12474 (2000).
  * [67] G. G. Batrouni, F. Hébert, and R. T. Scalettar, Phys. Rev. Lett. 97, 087209 (2006).
  * [68] B. Hetényi, B. Tanatar, and L. M. Martelo, Phys. Rev. B 93, 174518 (2016).
  * [69] G. Murthy, D. Arovas, and A. Auerbach, Phys. Rev. B 55, 3104 (1997).
  * [70] S. Wessel and M. Troyer, Phys. Rev. Lett. 95, 127205 (2005).
  * [71] D. Heidarian and K. Damle, Phys. Rev. Lett. 95, 127206 (2005); R. G. Melko, A. Paramekanti, A. A. Burkov, A. Vishwanath, D. N. Sheng, and L. Balents, Phys. Rev. Lett. 95, 127207 (2005); M. Boninsegni and N. Prokof’ev, Phys. Rev. Lett. 95, 237204 (2005).
  * [72] J. Gan, Y. Wen, and Y. Yu, Phys. Rev. B 75, 094501 (2007).
  * [73] F. Wang, F. Pollmann, and A. Vishwanath, Phys. Rev. Lett. 102, 017203 (2009).
  * [74] L. Pollet, J. D. Picon, H. P. Büchler, and M. Troyer, Phys. Rev. Lett. 104, 125302 (2010).
  * [75] X.-F. Zhang, S. Hu, A. Pelster, and S. Eggert, Phys. Rev. Lett. 117, 193201 (2016).
  * [76] A. Angelone, F. Mezzacapo, and G. Pupillo, Phys. Rev. Lett. 116, 135303 (2016).
  * [77] D. Yamamoto, I. Danshita, and C. A. R. Sá de Melo, Phys. Rev. A 85, 021601(R) (2012).
  * [78] M. Malakar, S. Ray, S. Sinha, and D. Angom, Phys. Rev. B 102, 184515 (2020).
  * [79] L. Santos, G. V. Shlyapnikov, and M. Lewenstein, Phys. Rev. Lett. 90, 250403 (2003).
  * [80] P. Nozières, J. Low Temp. Phys. 137, 45 (2004).
  * [81] J. Hertkorn, F. Böttcher, M. Guo, J. N. Schmidt, T. Langen, H. P. Büchler, and T. Pfau, Phys. Rev. Lett. 123, 193002 (2019).
  * [82] D. Pekker, B. Wunsch, T. Kitagawa, E. Manousakis, A. S. Sørensen, and E. Demler, Phys. Rev. B 86, 144527 (2012).
  * [83] D. S. Lühmann, Phys. Rev. A 87, 043619 (2013).
  * [84] J. Jin, A. Biella, O. Viyuela, L. Mazza, J. Keeling, R. Fazio, and D. Rossini, Phys. Rev. X 6, 031011 (2016).
  * [85] R. Bai, S. Bandyopadhyay, S. Pal, K. Suthar, and D. Angom, Phys. Rev. A 98, 023606 (2018).
  * [86] W. Krauth, M. Caffarel, and J. Bouchaud, Phys. Rev. B 45, 3137 (1992).
  * [87] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
  * [88] S. Ray, S. Sinha, and K. Sengupta, Phys. Rev. A 93, 033627 (2016).
  * [89] U. Pohl, S. Ray, and J. Kroha, Ann. Phys. (Berlin) 534, 2100581 (2022).
  * [90] K. W. Mahmud, E. N. Duchon, Y. Kato, N. Kawashima, R. T. Scalettar, and N. Trivedi, Phys. Rev. B 84, 054302 (2011).
  * [91] E. Altman and A. Auerbach, Phys. Rev. Lett. 89, 250404 (2002); S. D. Huber, E. Altman, H. P. Büchler, and G. Blatter, Phys. Rev. B 75, 085106 (2007).
  * [92] D. Pekker and C. M. Varma, Annu. Rev. Condens. Matter Phys. 6, 269 (2015).
  * [93] U. Bissbort, S. Götze, Y. Li, J. Heinze, J. S. Krauser, M. Weinberg, C. Becker, K. Sengstock, and W. Hofstetter, Phys. Rev. Lett. 106, 205303 (2011).
  * [94] M. Endres, T. Fukuhara, D. Pekker, M. Cheneau, P. Schauß, C. Gross, E. Demler, S. Kuhr, and I. Bloch, Nature (London) 487, 454 (2012).
  * [95] M. E. Fisher, M. N. Barber, and D. Jasnow, Phys. Rev. A 8, 1111 (1973).
  * [96] R. Roth and K. Burnett, Phys. Rev. A 68, 023604 (2003).
  * [97] V. L. Berezinskii, Sov. Phys. JETP 32, 493 (1971); J. M. Kosterlitz and D. J. Thouless, J. Phys. C: Solid State Phys. 6, 1181 (1973).
  * [98] F. Caleffi, M. Capone, C. Menotti, I. Carusotto, and A. Recati, Phys. Rev. Research 2, 033276 (2020).
  * [99] J. Zeiher, R. Van Bijnen, P. Schauß, S. Hild, J.-Y. Choi, T. Pohl, I. Bloch, and C. Gross, Nat. Phys. 12, 1095 (2016).
  * [100] P. Schauß, M. Cheneau, M. Endres, T. Fukuhara, S. Hild, A. Omran, T. Pohl, C. Gross, S. Kuhr, and I. Bloch, Nature (London) 491, 87 (2012).
  * [101] G. Schmid, S. Todo, M. Troyer, and A. Dorneich, Phys. Rev. Lett. 88, 167208 (2002).