Supersolid phases of Rydberg-excited bosons on a triangular lattice

TL;DR

This paper explores the phase diagram of Rydberg-excited bosons on a triangular lattice, revealing various supersolid phases and providing insights for experimental realization, with notable differences from square lattice behavior.

Contribution

It presents the first detailed numerical study of supersolid phases of Rydberg bosons on a triangular lattice, highlighting unique phase behaviors and proposing experimental parameters.

Findings

Identification of multiple supersolid phases on a triangular lattice.

Absence of checkerboard supersolid and superfluid-to-supersolid transition on the triangular lattice.

Comparison showing qualitative similarities and key differences with square lattice results.

Abstract

Recent experiments with ultracold Rydberg-excited atoms have shown that long-range interactions can give rise to spatially ordered structures. Observation of crystalline phases in a system with Rydberg atoms loaded into an optical lattice seems also within reach. Here we investigate a bosonic model on a triangular lattice suitable for description of such experiments. Numerical simulations based on bosonic dynamical mean-field theory reveal a rich phase diagram with different supersolid phases. Comparison with the results obtained for a square lattice geometry shows qualitatively similar results in a wide range of parameters, however, on a triangular lattice we do not observe the checkerboard supersolid. Moreover, unlike on a square lattice we did not find a phase transition from uniform superfluid to supersolid induced by increase of the hopping amplitude on a triangular lattice. Based on our results we propose an intuitive interpretation of the nature of different supersolid phases. We also propose parameters for the experimental realization.

Tables

Table 1. Table 1: Spanning vectors 𝐯 1 subscript 𝐯 1 \mathbf{v}_{1} and 𝐯 2 subscript 𝐯 2 \mathbf{v}_{2} and number of sites in the unit cell N u ​ c subscript 𝑁 𝑢 𝑐 N_{uc} of the structures considered in the B-DMFT calculations (except for the first one with N u ​ c = 3 subscript 𝑁 𝑢 𝑐 3 N_{uc}=3 , which was not considered explicitly but rather implicitly as a special case of the one with N u ​ c = 9 subscript 𝑁 𝑢 𝑐 9 N_{uc}=9 ). 𝐯 1 subscript 𝐯 1 \mathbf{v}_{1} and 𝐯 2 subscript 𝐯 2 \mathbf{v}_{2} are given in the basis of primitive vectors of a triangular lattice 𝐞 1 subscript 𝐞 1 \mathbf{e}_{1} and 𝐞 2 subscript 𝐞 2 \mathbf{e}_{2} , e.g. ( 3 , − 1 ) 𝐞 = 3 ​ 𝐞 1 − 𝐞 2 subscript 3 1 𝐞 3 subscript 𝐞 1 subscript 𝐞 2 (3,-1)_{\mathbf{e}}=3\mathbf{e}_{1}-\mathbf{e}_{2} . The primitive vectors in Euclidean space in units of the lattice spacing a 𝑎 a are 𝐞 1 = ( 1 , 0 ) subscript 𝐞 1 1 0 \mathbf{e}_{1}=(1,0) and 𝐞 2 = ( 1 2 , 3 2 ) subscript 𝐞 2 1 2 3 2 \mathbf{e}_{2}=(\frac{1}{2},\frac{\sqrt{3}}{2}) .

${𝐯}_1$	${(2,-1)}_{𝐞}$	${(2,0)}_{𝐞}$	${(3,-1)}_{𝐞}$	${(3,-1)}_{𝐞}$
${𝐯}_2$	${(1,1)}_{𝐞}$	${(0,2)}_{𝐞}$	${(0,2)}_{𝐞}$	${(1,2)}_{𝐞}$
$N_{uc}$	$3$	$4$	$6$	$7$
${𝐯}_1$	${(3,0)}_{𝐞}$	${(4,-2)}_{𝐞}$	${(4,-1)}_{𝐞}$	${(4,0)}_{𝐞}$
${𝐯}_2$	${(0,3)}_{𝐞}$	${(2,2)}_{𝐞}$	${(1,3)}_{𝐞}$	${(0,4)}_{𝐞}$
$N_{uc}$	$9$	$12$	$13$	$16$

Equations

$$\hat{H}={\hat{H}}_{kin}+{\hat{H}}_{vdW}+\sum _i\left({\hat{H}}_{loc,i}+{\hat{H}}_{R,i}\right).$$

$${\hat{H}}_{kin}=-J\sum _{⟨i,j⟩}\left({\hat{a}}_i^{†}{\hat{a}}_j^{}+\eta {\hat{b}}_i^{†}{\hat{b}}_j^{}\right).$$

$${\hat{H}}_{vdW}=\frac{V_{vdW}}{2}\sum _{i\ne j}\frac{{\hat{n}}_{e,i}{\hat{n}}_{e,j}}{|\mathbf{i}-\mathbf{j}{|}^6},$$

$$\begin{array}{rl}{\hat{H}}_{loc,i}& =\frac{U}{2}\left({\hat{a}}_i^{†}{\hat{a}}_i^{†}{\hat{a}}_i^{}{\hat{a}}_i^{}+2\lambda {\hat{a}}_i^{†}{\hat{b}}_i^{†}{\hat{b}}_i^{}{\hat{a}}_i^{}+\tilde{\lambda }{\hat{b}}_i^{†}{\hat{b}}_i^{†}{\hat{b}}_i^{}{\hat{b}}_i^{}\right)\\ & -\mu \left({\hat{n}}_{g,i}+{\hat{n}}_{e,i}\right),\end{array}$$

$${\hat{H}}_{R,i}=\frac{\Omega }{2}\left({\hat{b}}_i^{†}{\hat{a}}_i+{\hat{a}}_i^{†}{\hat{b}}_i\right)-\Delta {\hat{n}}_{e,i}.$$

$$\begin{array}{rl}{\hat{H}}_{vdW}& =\frac{V_{vdW}}{2}\sum _{i\ne j}\frac{{\hat{n}}_{e,i}{\hat{n}}_{e,j}}{|\mathbf{i}-\mathbf{j}{|}^6}\\ & \approx V_{vdW}\sum _{i\ne j}\left({\hat{n}}_{e,i}-\frac{⟨{\hat{n}}_{e,i}⟩}{2}\right)\frac{⟨{\hat{n}}_{e,j}⟩}{|\mathbf{i}-\mathbf{j}{|}^6}.\end{array}$$

$$z{J}_c=\frac{{\Omega }^2+{\Delta }^2-(2\mu +\Delta {)}^2}{4(\mu +\Delta )}.$$

$${\mathcal{G}}_i^{-1}(\mathrm{i}{\omega }_n)={\mathbf{G}}_i^{-1}(\mathrm{i}{\omega }_n)+{\mathbf{\Sigma }}_i(\mathrm{i}{\omega }_n).$$

$${\mathbf{\Psi }}_i=\left({\mathbf{g}}_i^0(0)-{\mathcal{G}}_i(0)\right){\mathbf{\Phi }}_i+\sum _{j:⟨i,j⟩}J{\mathbf{\Phi }}_j,$$

$${\mathcal{G}}_i(\mathrm{i}{\omega }_n)={\mathbf{g}}_i^0(\mathrm{i}{\omega }_n)-\sum _{j:⟨i,j⟩}{J}^2{\mathbf{G}}_j(\mathrm{i}{\omega }_n).$$

$${\tau }_c=\tau {\left(b\sum _i⟨{\hat{n}}_{e,i}⟩\right)}^{-1},$$

$${\Omega }_i=\Omega \exp \left[\frac{|\mathbf{i}-0{|}^2}{{\kappa }^2}\right],$$

$$\mathcal{V}=\frac{⟨{\hat{n}}_{g,i}⟩-⟨{\hat{n}}_{g,j}⟩}{⟨{\hat{n}}_{g,i}⟩+⟨{\hat{n}}_{g,j}⟩},$$

$$e_{\pm }=-\mu -\frac{\Delta \pm \sqrt{{\Delta }^2+{\Omega }^2}}{2}.$$

$$e_{+}\approx -\frac{1}{2}\left(1+\frac{{\Delta }_c}{\sqrt{{\Delta }_c^2+{\Omega }^2}}\right)\delta .$$

$$e_{vdW}\sim \rho {\int }_{r_c}^{\infty }\frac{r^{d-1}}{r^6}\mathrm{d}r=\rho \frac{r_c^{d-6}}{6-d}\sim {\rho }^{6/d}.$$

Full text

I Introduction

A supersolid is a phase with simultaneously broken $U(1)$ and translational symmetry of the system. Since the first time it was theoretically discussed Penrose and Onsager (1956); Andreev and Lifshitz (1969); Chester (1970); Leggett (1970) it has proven difficult to realize in experiment. So far it was only observed for ultracold bosons in optical cavities, where the light mode mediates long-range interaction Léonard et al. (2017); Landig et al. (2016); Klinder et al. (2015). Other paths to obtain supersolids are intensively studied. One promising experimental approach involves dipolar quantum gases loaded into an optical lattice Lahaye et al. (2009). The advantage of this approach is that the resulting system is highly tunable and accurately described by the extended Batrouni and Scalettar (2000); Góral et al. (2002) version of the Bose-Hubbard model Gersch and Knollman (1963); Fisher et al. (1989).

One of the first studies of the extended Bose-Hubbard model in the context of supersolids involved hard-core bosons on a triangular, frustrated lattice with nearest-neighbor interactions Wessel and Troyer (2005); Heidarian and Damle (2005); Melko et al. (2005); Boninsegni and Prokof’ev (2005). In these numerical quantum Monte Carlo (QMC) studies several different phases were observed, including a superfluid, an insulating density-wave, and two supersolid phases. Interestingly, the frustration of the lattice was found to be essential for the formation of the supersolid and for preventing phase separation, which is observed on the square lattice Batrouni and Scalettar (2000).

To obtain a supersolid phase on a square lattice with a hard-core constraint one needs long-range interaction van Otterlo et al. (1995); Ohgoe et al. (2011). Studies within mean-field Góral et al. (2002); Yi et al. (2007) and QMC Capogrosso-Sansone et al. (2010) have shown that supersolid phases exist in square and cubic lattice models with dipolar interactions between atoms. Long-range interaction also leads to other interesting phenomena, such as the appearance of multiple insulating density-wave phases with commensurate filling in the limit of small hopping amplitude Capogrosso-Sansone et al. (2010); Schachenmayer et al. (2010). The resulting phase diagram has features similar as the devil’s staircase in the Ising model Bak and Bruinsma (1982); Lauer et al. (2012); Rademaker et al. (2013).

Experimentally, long-range interactions of the van der Waals type can be realized, e.g. by exciting atoms to high principal quantum number Rydberg states Gallagher (1988); Singer et al. (2005). The advantage of this approach is that one can tune the strength of the long-range interaction through an appropriate choice of the Rydberg state. Rydberg atoms also have applications in other fields such as quantum information Saffman et al. (2010), understanding quantum critical behavior Weimer et al. (2008), molecule Bendkowsky et al. (2009); Niederprüm et al. (2015) and superatom Weber et al. (2015); Zeiher et al. (2015) formation, among many other Browaeys et al. (2016). Moreover, loading of Rydberg atoms into an optical lattice and the observation of emerging ordered structures have been achieved experimentally Schauß et al. (2012, 2015); Zeiher et al. (2016, 2017); Schauss (2018), although so far only in the frozen-limit of a deep lattice potential. In contrast, theoretical studies of the corresponding models have been performed both in the frozen-limit and for itinerant atoms. Main results include the study of self-organization of Rydberg excitations in a lattice Rademaker et al. (2013); Vermersch et al. (2015), phase diagrams and effects of finite hopping amplitude Lauer et al. (2012); Saha et al. (2014); Geißler et al. (2017); Li et al. (2018), spectral properties of different phases Geißler et al. (2018), and effects of dissipation Ray et al. (2016); Barbier et al. (2018). However, these studies focused on the square lattice geometry without considering effects of frustration.

In this work we aim at extending the study to the triangular lattice, frustrated with respect to formation of a checkerboard solid, and at understanding the influence of system’s geometry on the ground state by comparing our new results to those obtained for the square lattice Geißler et al. (2017). In Sec. II we introduce the model and briefly discuss our variant of the real-space bosonic dynamical mean-field theory (B-DMFT) method used in the simulations. In Sec. III we present and discuss the results. Sec. III.1 is aimed at understanding effects of the lattice geometry by studying the phase diagrams and properties of observed phases. In Sec. III.2 we compare results of B-DMFT and static mean-field theory to estimate the relevance of local quantum fluctuations. In Sec. III.3 we propose an experimental scheme for minimizing the destructive influence of dissipation induced by coupling to the environment. In Sec. IV our findings are summarized.

II System and method

II.1 Model

We choose a model suitable to describe experiments with bosonic alkali atoms, e.g. 87Rb, loaded into a triangular optical lattice Becker et al. (2010), and coupled (by an additional laser field) to a Rydberg state with high principal quantum number Schauß et al. (2012); Zeiher et al. (2016). For each lattice site $i$ we introduce bosonic annihilation operators $\hat{a}^{\phantom{{\dagger}}}_{i}$ of an atom in its ground state and $\hat{b}^{\phantom{{\dagger}}}_{i}$ of an atom in its highly excited Rydberg state. The corresponding Hamiltonian reads Saha et al. (2014); Geißler et al. (2017)

[TABLE]

The summation runs over the $N$ lattice sites of the system. In the end we take the thermodynamic limit of $N\to\infty$ assuming the system to be composed of periodically recurring unit cells of finite size $N_{uc}$.

$\hat{H}_{kin}$ represents the kinetic energy of atoms tunneling between neighboring lattice sites

[TABLE]

Here $J$ is the hopping amplitude, and $\eta$ represents the ratio between the hopping amplitude of excited state atoms to that of ground state atoms. $\langle i,j\rangle$ indicates summation over nearest neighbors $i$ and $j$. Is is useful to introduce the connectivity $z$ of the lattice, which is the number of nearest-neighbors for any site. For the triangular lattice $z=6$.

$\hat{H}_{vdW}$ represents the van der Waals interaction between two excited state atoms and is given by

[TABLE]

where $|\mathbf{i}-\mathbf{j}|$ is the Euclidean distance between lattice sites $i$ and $j$ divided by lattice spacing $a$, $\hat{n}_{e,i}=\hat{b}^{\dagger}_{i}\hat{b}^{\phantom{{\dagger}}}_{i}$ is the number operator at site $i$ for the excited bosons, $V_{vdW}$ is the van der Waals interaction strength, which is given by $V_{vdW}=C_{6}/a^{6}$ with $C_{6}$ being van der Waals coefficient Singer et al. (2005).

$\hat{H}_{loc,i}$ is a local part of the Hamiltonian (for site $i$) describing the chemical potential and the onsite interaction. It is given by

[TABLE]

with $\hat{n}_{g,i}=\hat{a}_{i}^{\dagger}\hat{a}^{\phantom{{\dagger}}}_{i}$. The parameters $U$, $\lambda U$ and $\tilde{\lambda}U$ describe the local interaction strength between two ground state atoms, ground state atom and excited state atom, and two excited state atoms, respectively. $\mu$ is the chemical potential of an external thermal reservoir, since we work in the grand canonical ensemble.

The last term in the Hamiltonian, the Rabi term $\hat{H}_{R,i}$, describes coupling between ground and excited state atoms, induced by the driving with an additional laser field. Within the rotating wave approximation (RWA) this contribution to the Hamiltonian is given by

[TABLE]

Here $\Omega$ is the Rabi frequency, and $\Delta$ the detuning of laser frequency from that of the atomic transition which we consider.

In the following we set $\hbar=k_{B}=1$ and use the Rabi frequency $\Omega$ as the unit of energy, unless stated otherwise. We assume the system is in thermal equilibrium at zero temperature.

In our model we set $\lambda,\tilde{\lambda}\gg 1$ leading to a hard-core constraint for excited state atoms Geißler et al. (2017). Rydberg atoms are susceptible to formation of molecules Manthey et al. (2015), which are not trapped by the lattice potential and therefore lead to a high two-body loss rate. This in turn leads to a hard-core constraint due to the quantum Zeno effect Misra and Sudarshan (1977); García-Ripoll et al. (2009); Vidanović et al. (2014).

We also set the value of $\eta=0$ which translates to immobile Rydberg atoms. The effect of non-vanishing $\eta$ was considered in Geißler et al. (2017) and only small changes in the values of observables were observed. This is not surprising as the excited atoms interact strongly via the van der Waals interaction, which leads either to a very small fraction of excited atoms or to crystalline order where kinetic processes are suppressed.

II.2 Method

We perform calculations with two methods: (i) a Gutzwiller (static) mean-field approximation, described in detail in Fisher et al. (1989); Barbier et al. (2018), and (ii) the bosonic dynamical mean-field theory (B-DMFT) Byczuk and Vollhardt (2008). Both methods are based on self-consistency and on mapping of the lattice problem onto a set of local impurity problems. In order to be able to do the latter we treat the non-local interaction term within the Hartree approximation Negele and Orland (1998)

[TABLE]

Moreover, both methods are implemented within a real-space approach, which allows to study arbitrary periodically recurring ordered structures.

Below we will outline the main steps of the B-DMFT approach, referring the reader to Geißler et al. (2017) for a more detailed discussion. The Gutzwiller mean-field technique may be viewed as a limiting case of B-DMFT and, therefore, it follows similar steps.

II.2.1 Frozen-limit

To efficiently perform the B-DMFT calculations we first need to predict what kind of self-organized structures may emerge in the system due to the long-range interaction. We therefore first perform calculations for the frozen gas with $J=0$. In this limit at unit filling one can map the problem onto an effective spin model Weber et al. (2015); Zeiher et al. (2016), which however is still not trivial to solve on an infinite two-dimensional triangular lattice. We therefore perform another simplification, assuming a negative value of the chemical potential $\mu<0$, which in the frozen-limit leads to a dilute crystal.

Owing to the negative chemical potential and zero temperature the bosons can reside in the lattice only when their energy is sufficiently lowered by the Rabi term $\hat{H}_{R}$. When this is the case, the ground state of the system will be a spatially periodic structure with optimal balance between the distribution of bosons in the system and the strength of interaction between them. One can efficiently find free energies of many metastable, spatially periodic states. Each such state is characterized by two spanning vectors $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$, see Fig. 1 (top). The bosons in the lattice reside only on the sites related by translations defined be these two vectors, forming a sublattice of the underlying triangular lattice. At each occupied lattice site there is exactly one boson $\langle\hat{n}_{g}+\hat{n}_{e}\rangle=1$ in a superposition of a ground and excited state. Comparing the free energies of these metastable ordered states, one can determine the ground state of the system. Structures associated in the frozen-limit with a ground state for a certain value of the detuning are considered in the B-DMFT calculation later on.

We note that in the frozen-limit a simple expression was found for the critical value of the detuning at which the system undergoes a phase transition to vacuum Geißler et al. (2017). This expression can be easily extended beyond the frozen-limit as the transition between the vacuum and a very dilute gas, in which van der Waals type interactions are negligible, can be treated as a single-particle problem (see also App. A). One finds a critical value of the hopping amplitude as a function of the chemical potential, detuning, Rabi frequency and connectivity

[TABLE]

We emphasize that the frozen-limit approach taken here is used primarily to predict most relevant structures for the further B-DMFT and static mean-field calculations. While in the frozen limit we do neglect certain orderings reported for lattice gas models that could not be described with just two Bravais vectors Lee et al. (2002); Rademaker et al. (2013), we still can recover some of them within B-DMFT (up to certain wavelengths of the structure) because there each site of the unit-cell is treated independently. We are not able to describe disordered, e.g. glassy, phases.

II.2.2 B-DMFT

The results of the frozen-limit allow to select the relevant ordered structures and thus to reduce the number of B-DMFT calculations by selecting only those pairs of spanning vectors $(\mathbf{v}_{1},\mathbf{v}_{2})$ which correspond to some ground state of the system in the frozen-limit. Each pair $(\mathbf{v}_{1},\mathbf{v}_{2})$ defines a unit cell with $N_{uc}$ sites that recurs periodically in the system. Within this unit cell a separate quantum impurity problem corresponds to each site. These impurity problems might have different parameters and solutions, resulting in different values of local observables, such as the condensate order parameter for ground $\phi_{i,g}=\langle\hat{a}_{i}\rangle$ and excited state $\phi_{i,e}=\langle\hat{b}_{i}\rangle$ bosons, expectation value of the occupation of ground $\langle\hat{n}_{g,i}\rangle$ and excited state $\langle\hat{n}_{e,i}\rangle$ bosons, connected local Green functions, self-energies, etc. As impurity solver within the B-DMFT calculations we apply the exact diagonalization method Hubener et al. (2009); Hu and Tong (2009); Geißler et al. (2017).

Within B-DMFT one needs to define a set of self-consistency equations Georges et al. (1996); Byczuk and Vollhardt (2008). The first one is given by the local Dyson equation, relating the local interacting connected Green function $\mathbf{G}$, local Weiss field $\mathbf{\mathcal{G}}$ and the self-energy $\mathbf{\Sigma}$. It reads

[TABLE]

Note that each object here is a $4\times 4$ matrix, since there are two components due to the Nambu notation for bosonic Green functions Byczuk and Vollhardt (2008) and two components due to two types of bosons (ground and excited state) in the lattice.

The second self-consistency equation involves the condensate order parameter and reads

[TABLE]

where $\mathbf{g}^{0}_{i}$ is the Green function of the non-interacting lattice site $i$ decoupled from the rest of the lattice, $\mathbf{\Psi}_{i}$ is the vector determining the condensate mean-field to which the impurity $i$ is coupled and $\mathbf{\Phi}_{j}$ is the vector determining the order parameter at site $j$ (which is calculated in the impurity problem). The summation runs over all $j$ which are nearest-neighbors of site $i$.

In standard B-DMFT the last self-consistency equation would be given by the lattice Dyson equation Byczuk and Vollhardt (2008). Here, however, due to the complexity of the problem and the large spatial structures considered we used a simpler one. In our approximate approach we determine the Weiss field according to

[TABLE]

Such self-consistency equation would become exact in the limit of the infinite connectivity Bethe tree Byczuk and Vollhardt (2008). In finite spatial dimension it amounts to neglecting: (i) the effect which removing a site from the lattice has on the lattice Green functions, (ii) correlations between different neighbors of the impurity. A similar (though not identical) self-consistency equation has been successfully used for lattices in finite dimensions in the context of real-time dynamics Strand et al. (2015). In App. B we elaborate on the effect of this approximation.

Having obtained the self-consistent solution one can use the local quantities calculated in the impurity problems to determine values of non-local quantities. Most importantly one can calculate the free energy Geißler et al. (2017). Note that we have included the chemical potential into the Hamiltonian (1) and are working at zero temperature, therefore the free energy per lattice site is given by $f=\langle\hat{H}\rangle/N$.

III Results

In this section we present the results of our calculations. In Sec. III.1 we set the system parameters to be comparable to those used in Geißler et al. (2017), where the same model on a square lattice has been studied with B-DMFT. This allows us to investigate how the triangular lattice geometry affects the behavior of the system. We also study the nature of different phases observed. In Sec. III.2 we compare the static mean-field and the B-DMFT results to estimate the significance of local quantum fluctuations. These two sections are aimed at giving a better understanding of phases emerging due to the competition of long-range interaction and kinetic processes on a triangular lattice. In Sec. III.3 we study a system with experimentally more feasible parameters. We investigate the possibility of observing supersolid phases in a triangular optical lattice with Rydberg atoms. To minimize dissipative effects we follow the idea suggested in Geißler et al. (2018) of using an inhomogeneous profile of the Rabi laser.

III.1 Phase diagram

We choose the following parameters of our system: $V_{vdW}=100\Omega$, $U=0.1\Omega$, $\mu=-0.25\Omega$. As discussed earlier we also set $\lambda=\tilde{\lambda}=10^{6}\gg 1$ and $\eta=0$. These parameters are the same as in Geißler et al. (2017), allowing for comparison of square and triangular lattices. Values of the hopping amplitude and the detuning are varied.

We first investigate the frozen-limit case $J=0$. In Fig. 1 (bottom) we show how the size of the unit cell of the ground state, given by its number of lattice sites $N_{uc}$, depends on the detuning. Below a critical value of $\Delta_{c}=-0.75$ the system is empty. As we increase $\Delta$ above $\Delta_{c}$ we observe a series of phase transitions between insulating ordered (density-wave) phases, resembling the devil’s staircase observed in the Ising model Bak and Bruinsma (1982); Lauer et al. (2012); Rademaker et al. (2013). Each has a different translational symmetry and size of the unit cell. For values of $\Delta$ close to $\Delta_{c}$ the unit cell is large, resulting in a very dilute system. As the value of the detuning is increased the density of bosons also increases. These results are qualitatively similar to the ones obtained for a square lattice Geißler et al. (2017). The trend line close to $\Delta_{c}$ follows $N_{uc}\sim(\Delta-\Delta_{c})^{-\frac{1}{3}}$. The exponent is determined by the spatial dimensionality of the system divided by the exponent in the interaction potential, cf. App. A.

Next we proceed to the discussion of the finite hopping $J>0$ case within the B-DMFT calculations. Out of a large set of unit-cells considered in the frozen-limit we have selected only the few smallest, relevant for the vicinity of $\Delta=0$, due to the computational complexity of B-DMFT calculations. They are listed in Tab. 1. According to the frozen-limit results the structures that were left out become relevant only in the narrow region of detuning $-0.75<\Delta<-0.7$ 111For much higher detuning one also expects a structure with $N_{uc}=3$ which was not included in Tab. 1. However, as it is a substructure of the one with $N_{uc}=9$ it can still be observed, similarly as $S_{6}$ in Fig. 2.. For other values of $\Delta$ the structures from Tab. 1 should be sufficient.

The phase diagram obtained within B-DMFT is shown in Fig. 2, while the density-wave patterns observed in different phases are shown in Fig. 3. The phases labeled as $\mathrm{DW}_{7}$, $\mathrm{DW}^{\prime}_{7}$, $\mathrm{DW}_{9}$ and $\mathrm{DW}_{12}$ are insulating while $\mathrm{SS}_{3}$, $\mathrm{SS}_{4}$, $\mathrm{SS}_{7}$, $\mathrm{SS}^{\prime}_{7}$, $\mathrm{SS}_{9}$, $\mathrm{SS}_{12}$ are supersolid. Lower index indicates number of sites in the unit cell $N_{uc}$ of the structure, cf. Tab.1. In the limit of small hopping amplitude $J=0.001$ we recover the results of the frozen-limit, as expected. The observed structures $\mathrm{DW}_{7}$ and $\mathrm{DW}_{9}$ (and $\mathrm{DW}_{12}$, not shown in Fig. 3) follow the trend presented in Fig. 1. Only the sites of a sublattice defined by vectors $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$ are occupied. On each of its sites there is a single boson, which is in a superposition between ground and excited state. The remaining sites of the lattice are nearly empty.

Next we consider the effect of increasing hopping amplitude. For $\Delta<\Delta_{c}=-0.75$ increasing $J$ leads to a phase transition from the vacuum to a homogeneous superfluid phase. The phase boundary obtained with B-DMFT agrees well with the expression (7). We note that below $\Delta_{c}$ unlike for the square lattice geometry, where a checkerboard supersolid was found Geißler et al. (2017), the triangular system does not exhibit any supersolid phase. For $\Delta>-0.6$ small values of the hopping amplitude have only a minor influence on the insulating phases, resulting in small shifts of the phase boundaries with increasing $J$. Further increase of $J$ eventually leads to a spontaneous breaking of the $U(1)$ symmetry in the system and a transition into one of many supersolid phases.

We observe that at higher values of $\Delta$ the supersolid phases immediately above their insulating counterparts have similar order, cf. Fig. 3, $\mathrm{DW}_{7}-\mathrm{SS}_{7}$ and $\mathrm{DW}_{9}-\mathrm{SS}_{9}$. In these phases, namely $\mathrm{SS}_{7}$ and $\mathrm{SS}_{9}$, a significant (when compared to ground state population) number of excited state bosons is present. We observe that both the local ground state condensate order parameter $\langle\hat{a}_{i}\rangle$ as well as local fluctuations of the occupation $\delta n_{i}^{2}=\langle(\hat{n}_{g,i}+\hat{n}_{e,i})^{2}\rangle-\langle\hat{n}_{g,i}+\hat{n}_{e,i}\rangle^{2}$ (see App.B) are much smaller on the sites with excited state bosons than on the surrounding sites with small excited state occupation 222The condensate of the excited state bosons is small throughout the system, with total fraction an order of magnitude smaller than that of the ground state bosons. Interestingly, even close to the $\mathrm{DW}_{9}$ phase it shows ordering later observed in $\mathrm{SS}_{3}$ supersolid.. We interpret this phase as a supersolid which consists of (i) frozen bosons being in a superposition of ground and excited state (with significant fraction of both) residing in the periodic sublattice of the original lattice, (ii) delocalized, condensed bosons, predominantly of the ground state nature, flowing without friction in the remaining lattice sites. This interpretation is also consistent with the shape of the supersolid–density-wave phase boundary for larger detunings. For a given ordered structure the delocalized bosons can slightly reduce their energy by a small admixture of the excited state, thus increasing the supersolid regime in the phase diagram. This energy reduction is however diminished with increasing fraction of excited state bosons in the “frozen” sites, due to strong non-local interaction forces. Hence with increasing $\Delta$ the phase boundary is shifted to higher values of $J$ as the fraction of excited bosons increases.

The supersolid phase $\mathrm{SS}_{12}$ appears to be of the same nature as $\mathrm{SS}_{7}$ and $\mathrm{SS}_{9}$ but with a different spatial structure. The $\mathrm{DW}_{12}$ phase (not depicted in Fig. 3), a $U(1)$ symmetric counterpart of the $\mathrm{SS}_{12}$ phase, appears for small $J$ around $\Delta\approx-0.5$, cf. Fig. 2. However, $\mathrm{SS}_{12}$ and its insulating counterpart are not immediately connected in the phase diagram due to the emergence of other phases discussed further in text.

Starting from the $\mathrm{SS}_{9}$ or $\mathrm{SS}_{12}$ phases and decreasing the detuning $\Delta$ or increasing the hopping amplitude $J$, we observe further phase transitions, cf. Fig. 2. Qualitatively different supersolids emerge, labeled as $\mathrm{SS}_{3}$ and $\mathrm{SS}_{4}$. These are characterized by the following features. The fraction of the excited state bosons is significantly (by approximately an order of magnitude) lower than in $\mathrm{SS}_{7}$, $\mathrm{SS}_{9}$ or $\mathrm{SS}_{12}$. The wavelength of the density-wave pattern is significantly smaller, with smaller distances between sites with non-vanishing fraction of excited state bosons. On these sites we have also observed an increase in condensate fraction and local fluctuations of the occupation, cf. App. B. Finally, we observe a larger density of atoms in the remaining, intermediate sites and more uniform distribution of the condensate order parameter. This behavior with decreasing values of $\Delta$ is opposite of what one would expect if one tried to apply here intuition gained from the frozen-limit. Because of these differences we conclude that this must be a qualitatively different type of supersolid, where we can no longer apply the interpretation of “frozen” sublattice sites occupied by the excited state bosons coexisting with condensed ground state bosons in between. These phases bear some resemblance to bubble supersolids Henkel et al. (2012); Cinti et al. (2014) in that the condensation originates form the sites with non-vanishing excited state fraction, in contrast to the $\mathrm{SS}_{7}$, $\mathrm{SS}_{9}$ and $\mathrm{SS}_{12}$ phases. We also suspect that the supersolids observed in $\mathrm{SS}_{3}$ and $\mathrm{SS}_{4}$ could be connected with the concept of defectons Andreev and Lifshitz (1969); Cinti et al. (2014), which is a condensation of defects (holes) in the ordered structure that are tunneling between different sites. This interpretation seems to be consistent with the observed features of the phases: (i) larger local fluctuations on the sublattice could be related to the presence and condensation of defectons, (ii) smaller wavelength of the observed pattern supports tunneling of defectons in opposition to the larger wavelength patterns, which are more favorable for the condensation of the ground state bosons on the intermediate sites (between the site of the sublattice). However, these arguments alone are not sufficient to confirm this interpretation unambiguously.

It is worth mentioning that there is a relation between $\mathrm{SS}_{9}$ and $\mathrm{SS}_{3}$ phases, as well as between $\mathrm{SS}_{12}$ and $\mathrm{SS}_{4}$ phases. In both cases the symmetry of the former phase can be viewed as a reduced version of the symmetry of the latter phase with $N_{uc}$ reduced by a factor of 3, cf. Fig. 3. However, the $\mathrm{SS}_{12}$ and $\mathrm{SS}_{4}$ phases, between which we observed a first order phase transition, separate the $\mathrm{SS}_{9}$ and $\mathrm{SS}_{3}$ phases from each other.

Regarding the two types of supersolids described above we note that similar observations were made in Li et al. (2018). There also two types were found, with one consisting of a supersolid of bare (ground state) species and crystalline phase of dressed (coupled to excited state) species. However, due to a different model the second type of supersolid in Li et al. (2018) is not of the same nature as observed here – one does not observe a significant reduction of the wavelength of the periodic structure but rather an increase.

We also observe a phase at intermediate values of the hopping amplitude $J\approx 0.038$ and in the vicinity of detuning $\Delta=0$ which is depicted in Fig. 3 and labeled $\mathrm{SS}^{\prime}_{7}$. It has similar properties to the $\mathrm{SS}_{3}$ and $\mathrm{SS}_{4}$ phases, which have a smaller wavelength than the phases at larger detuning, e.g. $\mathrm{SS}_{9}$, and for which condensation originates from the sites with non-vanishing excited state fraction, cf. App. B. The distinguishing features of $\mathrm{SS}^{\prime}_{7}$ are that the fraction of excited state bosons is significantly higher than in $\mathrm{SS}_{3}$ or $\mathrm{SS}_{4}$ and that we also observed an insulating counterpart of the $\mathrm{SS}^{\prime}_{7}$ phase, namely the $\mathrm{DW}^{\prime}_{7}$ phase. Both do not occur in the mean-field calculations, as discussed in the next section.

We note that the results presented here are qualitatively similar to the ones presented in Geißler et al. (2017) for a square lattice. A significant difference is visible here only near the boundary between superfluid and supersolid phases. On a square lattice this boundary separates superfluid from checkerboard supersolid, which extends to values of detuning below $\Delta_{c}$, where one can induce a phase transition from superfluid to supersolid by increasing (rather than decreasing) the hopping amplitude. On a triangular lattice the checkerboard supersolid cannot exist due to frustration and no supersolid phase exists below $\Delta_{c}$. Apart from this the two phase diagrams are similar. We attribute this to the fact that the spacing (wavelength) in the majority of structures observed here is larger than the lattice spacing. In this case the difference in geometry of the two lattices has a weaker impact. For these longer-wavelength structures it is actually the square lattice that becomes more frustrated with respect to the favored (due to van der Waals interaction) Wigner crystal formation Rademaker et al. (2013) than the triangular lattice.

III.2 Comparison with the static mean-field

In Fig. 4 we present the results of calculations performed within the static Gutzwiller mean-field approximation Fisher et al. (1989); Barbier et al. (2018). Upon comparison with the B-DMFT we notice that both methods give similar results. The main features of the phase diagrams agree well. Below we focus on the most relevant differences.

As the static mean-field approach favors ordered phases, we expect a phase transition between insulating and superfluid (supersolid) phases to apprear at lower values of the hopping amplitude. Indeed, comparing Fig. 2 with Fig. 4 we observe that the boundary of the insulating phases is shifted downwards. This effect is almost negligible for larger absolute values of the detuning and becomes relevant only in the region $|\Delta|\lessapprox 0.75$. This is also a region where an increasing number of phases compete in the system. It seems that only in this region local quantum fluctuations will significantly affect the system’s behavior.

Another discrepancy arises from the oversimplification of the insulating phases within the static mean-field approach. The boundaries between different phases with $U(1)$ symmetry do not depend on the hopping amplitude and their positions are uniquely defined by the detuning. In contrast, B-DMFT calculations show that finite hopping can induce a phase transition between two density-wave phases.

The above two observations are directly related to the most significant difference that we observed. Namely, in the static mean-field neither the supersolid $\mathrm{SS}^{\prime}_{7}$ nor the density-wave $\mathrm{DW}^{\prime}_{7}$ phase was observed, which should be present according to B-DMFT, cf. Fig. 2. This occurs in the region, where we observed largest discrepancies between the two methods in the values of the condensate order parameter and in properties of the insulating density-wave phase. We conclude that as the detuning $\Delta$ approaches critical value $\Delta_{c}=-0.75$ at intermediate values of the hopping amplitude, effects of local quantum fluctuations become significant (this is further backed up by comparing different self-consistency conditions, App. B). At this point we emphasize for clarity that non-local fluctuations are treated in both methods on the same level (Hartree mean-field). Therefore, we cannot make definite statements about their significance.

The last discrepancy between the results of B-DMFT and static mean-field can be observed in the extent of the $\mathrm{SS}_{4}$ and $\mathrm{SS}_{12}$ phases. In static mean-field it is slightly smaller than in B-DMFT due to larger extent of the $\mathrm{SS}_{9}$ phase. Nevertheless, the remaining features of the phase diagram are qualitatively accurately captured by the static mean field. As this method is significantly less demanding computationally, it is the best that we can do at this stage to get some insight into the critical region of $\Delta\approx\Delta_{c}$, where it is difficult to obtain converged B-DMFT results. We study this region of the phase diagram in Fig. 5.

As we approach $\Delta_{c}$ from above, we encounter a series of phase transitions forming a devil’s staircase (note the logarithmic scale), both of the insulating density-wave phases at small hopping amplitude as well as of the supersolid phases for larger hopping. When investigating the devil’s staircase in the supersolids we notice that increasing the hopping amplitude seems to favor longer wavelength structures shifting the devil’s staircase pattern to higher detunings. This is similar to what we observed for $\Delta\gtrapprox 1$, cf. Fig. 4. However, this trend seems to be reversed at the intermediate detuning of $\Delta\gtrapprox-0.68$ which coincides with the onset of $\mathrm{SS}_{3}$ phase (short wavelength) for larger hopping amplitude. We suspect that this might be a feature emerging due to the competition between the two general types of supersolid discussed in Sec. III.1 and could lead to a multicritical point around $\Delta\approx-0.68$ and $J\approx 0.022$. However, investigating this region within a more accurate B-DMFT method requires improvements of our implementation of the method and goes beyond the scope of this work.

III.3 Finite-size system with inhomogeneous Rabi frequency

In Sec. III.1 we have established a relation between the behavior of long-range interacting bosons on the triangular and square lattices. In order to make this comparison, we have chosen the same parameters as used in Geißler et al. (2017). However, these values of the parameters are not optimally suited for experimental realization of the model (1) with Rydberg atoms loaded into an optical lattice. Relatively large values of the hopping amplitude $J$ and small values of the van der Waals interaction $V_{vdW}$ with respect to the Rabi frequency would require using a Rydberg excited state with low principal quantum number $n\sim 16$. These states have a short lifetime due to spontaneous emission and dephasing processes induced by the black-body radiation Beterov et al. (2009). Below we propose more realistic parameters, which were chosen based on the observations made in the previous sections and for which a supersolid phase could be observed.

The main two sources of dissipation in the system are the spontaneous emission and black-body induced dephasing Beterov et al. (2009). The lifetime $\tau$ of an excited Rydberg state depends approximately as $\tau\sim n^{3}$ on the principal quantum number $n$. It is therefore advantageous to increase $n$.

Another issue one needs to overcome is the avalanche dephasing Zeiher et al. (2017); Boulier et al. (2017); Geißler et al. (2018). A black-body radiation-induced transition of a single atom to another Rydberg state might trigger an avalanche of atom loss from the system. The average time after which such a process occurs is given by Boulier et al. (2017)

[TABLE]

where $b$ is the branching ratio of the excited state. $\tau_{c}$ is inversely proportional to the total number of excited state atoms in the system. Therefore, the best candidate for experimental observation of a supersolid phase is the $\mathrm{SS}_{3}$ phase where the excited state fraction is very low. However, since in (11) the total number of Rydberg atoms in the system appears, rather than their density, we need to consider relatively small system sizes. This leads to the further issue of increased Rydberg fraction at sharp edges of the finite size system, such as shown in Vermersch et al. (2015); Geißler et al. (2018). In order to avoid this problem we additionally choose a Gaussian profile of the Rabi laser with a narrow waist on the order of several $\mu$m, as suggested in Geißler et al. (2018), given by

[TABLE]

with $\mathbf{0}$ corresponding to the position of the center of the system and $\Omega_{i}$ an effective Rabi frequency at site $i$.

Taking the above restrictions into account we consider a system with bosonic 87Rb atoms loaded into a two dimensional triangular optical lattice, e.g. such as described in Becker et al. (2010) with lattice spacing $a\approx 0.5\mu$m. For such a system one should be able to achieve a hopping amplitude of approximately $\sim 0.45\hbar$kHz Becker et al. (2010). The local interaction is a tunable parameter, which we set to $U=1\hbar$kHz. We choose to couple the ground state to an excited $|26S\rangle$ Rydberg state by the Rabi term. Using this $|26S\rangle$ state on a lattice with spacing $a\approx 0.5\mu$m gives a van der Waals interaction strength on the order of $V_{vdW}\approx 1600\hbar$MHz Singer et al. (2005) and mean lifetime of $\tau\approx 10\mu$s Beterov et al. (2009). Rabi frequency and detuning on the order of $\Omega\approx 0.1\hbar$MHz, $\Delta=-0.4\hbar$MHz should also be feasible experimentally. The remaining issue is to focus the Rabi laser such that it has a Gaussian profile (12) with $\kappa=3.5\mu$m.

Setting everything in units of $\Omega$ and $a$ we obtain the parameters of the simulation to be the following: $J=0.0045\Omega$, $\mu=-0.025\Omega$, $U=0.01\Omega$, $\Delta=-4\Omega$, $V_{vdW}=1.6\times 10^{4}\Omega$ and $\kappa=7a$. The resulting density pattern is shown in Fig. 6.

In the center of the system we observe the same structure as in the $\mathrm{SS}_{3}$ phase. Its visibility, defined as

[TABLE]

where $i$ corresponds to the site with maximal occupation and $j$ to the nearest-neighbor of $i$, has the value $\mathcal{V}\approx 0.175$.

Lastly we note that the average total number of atoms in the system is $\sum_{i}\langle\hat{n}_{g,i}+\hat{n}_{e,i}\rangle\approx 113$ while the average total number of Rydberg excitations is $\sum_{i}\langle\hat{n}_{e,i}\rangle\approx 0.0373$. Together with the branching ratio on the order of $b\approx 0.1$ Boulier et al. (2017) the average time after which the avalanche is set off is $\tau_{c}\approx 3$ms, which is comparable with the characteristic timescale of the hopping process $t\sim\hbar/J$, and thus is promising for experimental realization. The remaining challenges are to achieve low enough temperatures of the system and to determine whether dissipative effects other than avalanche dephasing can destroy the order in the supersolid phase. The latter has been investigated in Ref. Barbier et al. (2018) and it seems that even for relatively large dissipation strength the order prevails for times on the order of hundreds of $\mu$s if the fraction of excited state atoms is small, similarly as in the case of avalanche dephasing.

IV Conclusions

In this work we have studied the effect of frustration on the formation of crystalline and supersolid states in the extended Bose-Hubbard model with two bosonic species, one itinerant and one subject to two-body van der Waals long-range interaction. The two species are also coupled by a Rabi term and local interaction. We have focused on a system at zero temperature, without including explicitly the effects of coupling to the environment. We have also used the Hartree mean-field approximation to decouple the long-range interaction term.

Within B-DMFT the same model has been previously investigated on a two-dimensional square lattice Geißler et al. (2017). Here we solve the problem with two methods: B-DMFT and the Gutzwiller static mean-field approach. Comparison of the B-DMFT results for the two lattice geometries allowed us to determine the effect of the (frustrated) triangular geometry on ordered states. Comparison of B-DMFT results to those of the static mean-field approach allowed us to estimate the significance of local quantum fluctuations.

We have obtained a rich phase diagram, including: insulating density-wave, superfluid and supersolid phases. We observed that the phase diagram on the triangular lattice is qualitatively similar to the one observed for a system with square lattice geometry. Within the parameter regimes considered the only significant discrepancy is the absence of the checkerboard ordered supersolid and of the supersolid below a critical value of the detuning $\Delta_{c}$ (determined in the frozen-limit). The similarity of the results for the two geometries can be attributed to the low density of atoms and the large wavelength of the observed ordered structures, when compared to the lattice spacing.

Comparison of the B-DMFT and static mean-field methods shows an overall good agreement between the two approaches. The differences are limited to small regions of the phase diagram and small phase-boundary shifts.

We have also studied the model on a finite-size lattice with Gaussian profile of the Rabi term, where we have chosen experimentally convenient parameters. We have found a supersolid phase with a low fraction of excited atoms and a visible spatial modulation of the density. We believe this to be the most promising approach for experimental realization of supersolid phases with Rydberg atoms, minimizing the effects of dissipation.

Acknowledgements.

Support by the Deutsche Forschungsgemeinschaft via DFG SPP 1929 GiRyd and the high-performance computing center LOEWE-CSC is gratefully acknowledged. The authors also acknowledge useful discussions with C. Groß, S. Hollerith, W. Li, Y. Li, S. Whitlock and J. Zeiher.

Appendix A Scaling of $N_{uc}$ close to $\Delta_{c}$ in the frozen-limit

Below we give a simplified argumentation for the dependence of the size of the unit cell $N_{uc}$ on the detuning close to the critical detuning strength $\Delta_{c}$ in the frozen-limit $J=0$. We assume the system is $d$-dimensional. We first consider the energy gain due to adding a single particle into an empty system. In such case the only relevant energy scales in the Hamiltonian (1) are the chemical potential $\mu$, the detuning $\Delta$ and the Rabi frequency $\Omega$. In order to have a finite value of $\Delta_{c}$, which marks the transition of the system to vacuum, we set $\mu<0$. Finding the single-particle eigenstates in such case is simple and their energies are given by

[TABLE]

The low energy state is given by $e_{+}$, and when $e_{+}<0$ it is energetically favorable to put the particles in the system. Thus the condition $e_{+}=0$ determines the critical value of detuning $\Delta_{c}$ (e.g. in case of $\Omega=1$ and $\mu=-0.25$ we get $\Delta_{c}=-0.75$, cf. Sec. III.1). Note that this argumentation can be extended beyond the frozen-limit of $J=0$. By taking a completely delocalized single particle state of the ground state bososns, giving $-zJ$ contribution to the energy, one can obtain critical value of the hopping amplitude $J_{c}$ as a function of the chemical potential, the detuning and the Rabi frequency, yielding (7). We further consider the energy gain per particle in the vicinity of $\Delta_{c}$, with $\Delta=\Delta_{c}+\delta$ and $0<\delta\ll 1$. Expanding (14) up to first order in $\delta$ we get

[TABLE]

Therefore, energy gain per particle due to adding particles in the system is proportional to $\delta=(\Delta-\Delta_{c})$.

However, upon adding particles into the system we increase the potential energy due to the van der Waals interaction. We therefore need to estimate the energy cost due to having a certain density of particles in the system. We assume that the particles form a uniform Wigner crystal with certain density $\rho\sim 1/N_{uc}$. In such case the average distance between particles is given by $r_{c}\sim\rho^{-1/d}$. The energy per particle due to the van der Waals interaction can be estimated by integral

[TABLE]

If we now require the energy cost of interaction to be compensated by the energy gain due to the Rabi frequency we obtain $(\Delta-\Delta_{c})\sim\rho^{6/d}\sim N_{uc}^{-6/d}$, which for $d=2$ gives $N_{uc}\sim(\Delta-\Delta_{c})^{-1/3}$, cf. Fig. 1.

Appendix B Supplementary results

B.1 Self-consistency test

In order to determine the accuracy of our assumption regarding the self-consistency condition (10) described in Sec. II.2 we have implemented the full self-consistency condition and performed calculations for limited range of parameters. We have chosen to set the detuning to $\Delta=0$ and vary the hopping amplitude $J$ as this gives a cross section of the most interesting part of the phase diagram shown in Fig. 2. We compare the average occupation of the excited state bosons as a function of $J$. The result is shown in Fig. 7.

One can see that for the majority of values of the hopping amplitude $J$ the two self-consistency equations yield quantitatively comparable results. The major difference appears in the region of $J\in[0.031,0.03925]$. In the simplified self-consistency the extent of $\mathrm{DW}^{\prime}_{7}$ phase is much smaller while the extent of $\mathrm{SS}^{\prime}_{7}$ phase is slightly larger. In the latter phase the simplified approach yields also larger rate of change of the average occupation of excited state bosons with increasing $J$. We note that while the extent of these phases is affected, the general features of the phase diagram remain unaffected. We do not observe significant differences, e.g. appearance of new types of phases. Tests for other values of $\Delta$, not shown here, confirmed this conclusion.

B.2 Condensate order parameter and local fluctuations

In Fig. 8 we present additional results showing local condensate fraction $\phi_{i,g}=\langle\hat{a}_{i}\rangle$ and local fluctuations of site occupation $\delta n^{2}_{i}=\langle\hat{n}_{i}^{2}\rangle-\langle\hat{n}_{i}\rangle^{2}$, for different sites $i$. In the phases $\mathrm{SS}_{7}$, $\mathrm{SS}_{9}$ and $\mathrm{SS}_{12}$ we observe that both the condensate order parameter and local fluctuations are significantly suppressed at the sites occupied by the excited state bosons. This is consistent with our interpretation that in these phases we observe frozen bosons on selected sites with non-vanishing excited state fraction and condensed bosons in the intermediate sites, which are responsible for the superflow. A qualitatively different behavior is observed for the $\mathrm{SS}_{3}$, $\mathrm{SS}_{4}$ and $\mathrm{SS}^{\prime}_{7}$ phases. There the condensate fraction and local fluctuations are actually larger at the sites with a significant excited state bosons fraction. This property is reminiscent of the bubble supersolids observed in Ref. Henkel et al. (2012); Cinti et al. (2014), although here we work in the significantly different regime of small detuning. We note that the spatial modulation of the two observables considered here is much smaller in the $\mathrm{SS}_{3}$ and $\mathrm{SS}_{4}$ phase than in the $\mathrm{SS}^{\prime}_{7}$ phase.