Supersolidity of lattice Bosons immersed in strongly correlated Rydberg dressed atoms

TL;DR

This paper demonstrates that interspecies interactions in a two-component lattice boson system with Rydberg dressing can induce and stabilize supersolid phases, revealing new quantum states accessible in current experiments.

Contribution

It introduces a model of two-component lattice bosons with Rydberg-induced long-range interactions and shows how interspecies interactions stabilize supersolidity, a novel mechanism in this context.

Findings

Interspecies onsite interactions stabilize supersolid phases.

Two types of supersolids are identified: crystalline and immersed in quantum phases.

Interspecies interactions cause roton-like instabilities in the bare species.

Abstract

Recent experiments have illustrated that long range two-body interactions can be induced by laser coupling atoms to highly excited Rydberg states. Stimulated by this achievement, we study supersolidity of lattice bosons in an experimentally relevant situation. In our setup, we consider two-component atoms on a square lattice, where one species is weakly dressed to an electronically high-lying (Rydberg) state, generating a tunable, soft-core shape long-range interaction. Interactions between atoms of the second species and between the two species are characterized by local inter- and intra-species interactions. Using a dynamical mean-field calculation, we find that interspecies onsite interactions can stabilize a pronounced region of supersolid phases. This is characterized by two distinctive types of supersolids, where the bare species forms supersolid phases that are immersed in strongly correlated quantum phases, i.e. a crystalline solid or supersolid of the dressed atoms. We show that the interspecies interaction leads to a roton-like instability in the bare species and therefore is crucially important to the supersolid formation. We provide a detailed calculation of the interaction potential to show how our results can be explored under current experimental conditions.

Figures

Equations

$$\hat{H}=$$

$${\hat{H}}_{\mathrm{e}}=$$

$$S_{\text{imp}}^{(i)}$$

$$\displaystyle\bm{\mathcal{G}}^{-1}_{0,\sigma\sigma^{\prime}}(\tau-\tau^{\prime})\equiv-\left(\begin{array}[]{cc}(\partial_{\tau^{\prime}}-\mu_{\sigma})\delta_{\sigma\sigma^{\prime}}+t^{2}\sum\limits_{\langle 0i\rangle,\langle 0j\rangle}G_{\sigma\sigma^{\prime},ij}^{1}(\tau,\tau^{\prime})&t^{2}\sum\limits_{\langle 0i\rangle,\langle 0j\rangle}G^{2}_{\sigma\sigma^{\prime},ij}(\tau,\tau^{\prime})\\ t^{2}\sum\limits_{\langle 0i\rangle,\langle 0j\rangle}{G^{2}_{\sigma\sigma^{\prime},ij}}^{*}(\tau^{\prime},\tau)&(-\partial_{\tau^{\prime}}-\mu_{\sigma})\delta_{\sigma\sigma^{\prime}}+t^{2}\sum\limits_{\langle 0i\rangle,\langle 0j\rangle}G_{\sigma\sigma^{\prime},ij}^{1}(\tau^{\prime},\tau)\end{array}\right),$$

$${\varphi }_{i,\sigma }(\tau )\equiv ⟨b_{i,\sigma }(\tau ){⟩}_0$$

$$G_{\sigma {\sigma }^{\prime },ij}^1(\tau ,{\tau }^{\prime })$$

$$G_{\sigma {\sigma }^{\prime },ij}^2(\tau ,{\tau }^{\prime })$$

$${\hat{H}}_A=$$

$$\frac{1}{2}\sum _{i\ne j}{V}_{ij}{\hat{n}}_{i,d}{\hat{n}}_{j,d}\approx \sum _{i\ne j}{V}_{ij}⟨{\hat{n}}_{i,d}⟩({\hat{n}}_{j,d}-\frac{1}{2}⟨{\hat{n}}_{i,d}⟩)$$

$$G_{\mathrm{imp},\sigma {\sigma }^{\prime }}^1(i{\omega }_n)$$

$$G_{\mathrm{imp},\sigma {\sigma }^{\prime }}^2(i{\omega }_n)$$

$${\mathbf{\Delta }}_{\sigma {\sigma }^{\prime }}(i{\omega }_n)$$

$${\Delta }_{\sigma {\sigma }^{\prime }}^1(i{\omega }_n)$$

$${\Delta }_{\sigma {\sigma }^{\prime }}^2(i{\omega }_n)$$

$${\mathbf{\Sigma }}_{\mathrm{imp},\sigma {\sigma }^{\prime }}(i{\omega }_n)=(i{\omega }_n{\sigma }_z+{\mu }_{\sigma }){\delta }_{\sigma {\sigma }^{\prime }}+{\mathbf{\Delta }}_{\sigma {\sigma }^{\prime }}(i{\omega }_n)-{\mathbf{G}}_{\mathrm{imp},\sigma {\sigma }^{\prime }}^{-1}(i{\omega }_n).$$

$$\mathbf{G}(i{\omega }_n{)}^{-1}={\mathbf{G}}_0(i{\omega }_n{)}^{-1}-\mathbf{\Sigma }(i{\omega }_n).$$

$${\mathbf{G}}_0(i{\omega }_n{)}^{-1}=(\mu +i{\omega }_n)1-\mathbf{t}.$$

$$\Big{(}\bm{\mathcal{G}}^{(i)}_{0,\sigma\sigma^{\prime}}(i\omega_{n})\Big{)}^{-1}=\Big{(}\mathbf{G}^{(i)}_{\sigma\sigma^{\prime}}(i\omega_{n})\Big{)}^{-1}+\mathbf{\Sigma}^{(i)}_{\sigma\sigma^{\prime}}(i\omega_{n}),$$

$${\hat{H}}_{\mathrm{kin}}=-\sum _{⟨ij⟩,\sigma }{t}_{\sigma }(b_{i,\sigma }^{†}{b}_{j,\sigma }+H.c.).$$

$$E_{\mathrm{kin}}$$

$$=\sum _{ij,\sigma }{t}_{ij}^{\sigma }\left(\underset{ϵ\to 0^{+}}{\lim }\sum _{n=-\infty }^{\infty }\frac{e^{i{\omega }_nϵ}}{\beta }{G}_{ji,\sigma }(i{\omega }_n)-{\varphi }_{i,\sigma }^{\ast }{\varphi }_{j,\sigma }\right)$$

$${\mathbf{G}}_i^C(i{\omega }_n{)}^{-1}$$

$$[{\mathbf{G}}_R^C(i{\omega }_n{)}^{-1}{]}_{ij}$$

$$\sum _j[{\mathbf{G}}_R^C(i{\omega }_n{)}^{-1}{]}_{ij}[\mathbf{G}(i{\omega }_n){]}_{ji}=\sum _j\left[t_{ij}{1}_2-{\delta }_{ij}\left({\mathbf{\Delta }}_i(i{\omega }_n)-{\mathbf{G}}_i(i{\omega }_n{)}^{-1}\right)\right][\mathbf{G}(i{\omega }_n){]}_{ji}$$

$$\sum _j{t}_{ij}[\mathbf{G}(i{\omega }_n){]}_{ji}={\mathbf{\Delta }}_i(i{\omega }_n){\mathbf{G}}_i(i{\omega }_n).$$

$$E_{\mathrm{kin}}=$$

$$-\sum _{ij\sigma }{t}_{ij}{\varphi }_{i,\sigma }^{\ast }{\varphi }_{j,\sigma }-\frac{\text{Tr}\left[{\mathbf{\Delta }}_{\sigma ,i}(0){\mathbf{G}}_{\sigma ,i}(0)\right]}{2\beta }.$$

$$E=E_{\mathrm{kin}}+E_{\mathrm{int}}.$$

$$E_{\mathrm{int}}=\frac{1}{2}\sum _{i,\sigma {\sigma }^{\prime }}{U}_{\sigma {\sigma }^{\prime }}{\hat{n}}_{i,\sigma }({\hat{n}}_{i,{\sigma }^{\prime }}-{\delta }_{\sigma {\sigma }^{\prime }})+\sum _{i<j}{V}_{ij}{\hat{n}}_{i,d}{\hat{n}}_{j,d}.$$

$${\hat{H}}_{\mathrm{e}}=$$

$$\tilde{H}=-\sum _{\vec{k}}[\mu +2t(\cos k_{x}a+\cos k_{y}a)]b_{\vec{k}}^{†}{b}_{\vec{k}}+\frac{U}{2N}\sum _{{\vec{k}}_1{\vec{k}}_2{\vec{k}}_3}{b}_{{\vec{k}}_1}^{†}{b}_{{\vec{k}}_3}^{†}{b}_{{\vec{k}}_3+{\vec{k}}_2}{b}_{{\vec{k}}_1-{\vec{k}}_2}+U_1\sum _{\{\vec{k}\}}{b}_{\vec{k}}^{†}{b}_{\vec{k}},$$

$$\tilde{H}\approx E_0-\sum _{\vec{k}\ne 0}\left[\mu +2t(\cos k_{x}a+\cos k_{y}a)-2U{\bar{n}}_b\right]b_{\vec{k}}^{†}{b}_{\vec{k}}+\frac{U{\bar{n}}_b}{2}\sum _{\vec{k}\ne 0}(b_{\vec{k}}{b}_{-\vec{k}}+b_{-\vec{k}}^{†}{b}_{\vec{k}}^{†})+U_1\sum _{\{\vec{k}\ne 0\}}{b}_{\vec{k}}^{†}{b}_{\vec{k}},$$

$$\tilde{H}\approx E_0+\sum _{\vec{k}\ne 0}\left[{\epsilon }_k+U{\bar{n}}_b\right]b_{\vec{k}}^{†}{b}_{\vec{k}}+\frac{U{\bar{n}}_b}{2}\sum _{\vec{k}\ne 0}(b_{\vec{k}}{b}_{-\vec{k}}+b_{-\vec{k}}^{†}{b}_{\vec{k}}^{†}),$$

Full text

A supersolid is a translational symmetry breaking superfluid occurring in a solid. It was predicted to exist in bulk helium over forty years ago supersolid , but its observation has remained a challenge Kim . To reach supersolidity, one typically relies on long-range two-body interactions to break the translational invariance of a homogeneous system. Recent experiments have observed supersolid orders where translational symmetry is broken by cavity photon assisted cavity or spin-orbit coupling enabled SOC momentum transfer.

To achieve supersolids induced purely by two-body interactions, enormous efforts have been spent on polar molecules Ye ; Ye_theory , magnetic Cy_Dr_Ey and Rydberg atoms single_component ; single_eBH , due to the available long-range atom-atom interaction as well as high precision control over their internal and motional states. However, a current challenge is that theoretical proposals typically examine regimes that are difficult to achieve experimentally.

In this work, we study supersolids of a two-species bosonic mixture on a two-dimensional (2D) square lattice, where one of the species is weakly coupled to an electronically high-lying (Rydberg) state by an off-resonant laser (the level scheme is depicted in Fig. 1a). Uniquely, this setting is recently realized experimentally at Munich Rydberg in the study of Rydberg dressed spin dynamics spin_lattice . The coupling laser induces strong and long-range interactions between Rydberg dressed atoms on distances well beyond typical lattice-site spacings (see Fig. 1b), whose strength and sign can be controlled by the laser (i.e. detuning and Rabi frequencies) and the choices of Rydberg states dressed . The resulting Bose-Hubbard model features a long-range interaction between dressed atoms while interactions between atoms of the two different species and of the bare species are short ranged.

Employing real-space bosonic dynamical mean-field theory (RBDMFT), we find that the system undergoes a series of many-body phases, including Mott insulator (MI), ordered density wave (DW), supersolid (SS) and superfluid (SF) phases. A key result is that the interspecies interaction enables supersolid phases of the bare species in regions where the dressed atoms are in DW or SS phases (an example for a DW is depicted in Fig. 1c). Using Bogoliubov theory, we reveal that a roton-like instability emerges due to the interspecies interaction (see Fig. 1d), which signifies a SF to SS transition SF_SS . Our results open a new route to enhance the formation of SS phases through the Rydberg dressing in two-component atomic gases.

The Hamiltonian– In sufficiently deep lattices, our setting is described by a single band, two-component Bose-Hubbard model,

[TABLE]

where the single site Hamiltonian $\hat{H}_{i}=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}U_{\sigma\sigma^{\prime}}\hat{n}_{i\sigma}(\hat{n}_{i\sigma^{\prime}}-\delta_{\sigma\sigma^{\prime}})-\sum_{\sigma}\mu_{\sigma}\hat{n}_{i\sigma}$. $\langle i,j\rangle$ represents the nearest neighbour sites $i,j$. Index $\sigma(\sigma^{\prime})=b,d$ denotes bare, and dressed states, respectively. $\hat{b}^{\dagger}_{i\sigma}$ ($\hat{b}_{i\sigma}$) and $\hat{n}_{i\nu}=\hat{b}^{\dagger}_{i\nu}\hat{b}_{i\nu}$ are the bosonic creation (annihilation) operator for species $\sigma$ and atomic density at site $i$. $t$ and $\mu_{\sigma}$ determine the hopping rate and chemical potential for the two bosonic species. We assume the hopping rates are identical for both species identical_hopping . $U_{\sigma\sigma^{\prime}}$ denotes the inter- and intra-species short-range (onsite) interactions, which can be tuned via e.g. Feshbach resonances Feshbach or state-dependent optical lattices state_lattice . The long-range interaction between site $i$ and $j$ is $V_{ij}\equiv V/[(a/R_{c})^{6}(i-j)^{6}+1]$, where $V=\tilde{C}_{6}/R_{c}^{6}$ characterises the long-range interaction at a distance $R_{c}$. $\tilde{C}_{6}$, $R_{c}$ and $a$ are the effective dispersion coefficient, soft-core radius, and lattice constant, respectively. In the following, we choose the intraspecies short-range interaction $U_{b,d}\equiv U$, which also sets the unit of energy. Details of these parameters will be given towards the end of the paper.

To determine the ground state phases, we use RBDMFT to capture both higher order quantum fluctuations, strong correlations and arbitrary long-range order in a unified framework RBDMFT ; BDMFT_2009 . It provides a nonperturbative description of many-body systems in two and three spatial dimensions (the method is discussed in the supplementary material.). In the calculations, we typically consider the lattice size as large as $N_{\rm lat}=48\times 48$ sites and an experimentally relevant soft-core radius $R_{c}=3a$ dressed . The superfluidity is characterised by the condensate order parameter $\phi_{\sigma}\equiv\langle\hat{b}_{\sigma}\rangle$, and crystalline order by the real-space density distribution $n_{i\sigma}=\langle\hat{n}_{i\sigma}\rangle$ and total density $n_{i}\equiv n_{ib}+n_{id}$. The coexistence of both condensate and crystalline order parameters gives the supersolid phase. Note that a similar model using dipolar gases has been numerically investigated using a mean-field Gutzwiller approach and by considering only the nearest-neighbor part of the dipolar interactions NN . In our calculations, we take into account the whole range of the interaction potential (see appendix for a comparison of the two systems).

Many-body ground state phase diagram— The main results are summarized in the phase diagram shown in Fig. 2. Depending on the parameters, the two-component system can have five different phases, i.e. Mott insulator (MI), ordered density wave (DW), two types of supersolid (SS1 and SS2), and superfluid (SF). In the following, we will discuss features of these phases for unit filling $n_{jd}+n_{jb}=1$ (see appendix for results at other fillings).

We start with the so-called strong coupling limit when $U_{\sigma\sigma^{\prime}}\gg t$, where the 2D system favours MI phases with uniform total particle densities. Crystalline orders in the MI region can be changed by varying the two-body interactions (i.e. $V/U$). One example is depicted in Fig. 2a, which shows relative densities and crystalline structures. Furthermore, when one increases $V/U$ continuously, the filling fractions $f_{d}\equiv\sum_{i}n_{id}/N_{\rm lat}$ of the dressed species can form a devil’s staircase structure (Fig. 3a). An open question here is whether the staircase in this 2D system is complete. In 1D lattice systems, the devil’s staircase and its completeness one_dimension have been extensively studied staircase . Moreover, there are very small regions occupied by DW phases (with a non-uniform total density). Due to that, the corresponding discussion will be given in the supplementary material.

When the hopping rate increases, we observe a pronounced region of supersolids. The bare state first enters the supersolid phase (SS1) from an insulating phase, while the dressed species is still crystallized in this case (one example is depicted in Fig. 2b). Further increasing $t$, both species are in supersolid phases (SS2), as shown in Fig. 2c, where non-zero peaks appear for both species in addition to zero-momentum condensate, indicating the coexistence of non-trivial diagonal long-range order and off-diagonal long-range order associated with phase coherence. A large supersolid region indicates a higher chance for directly observing these phases in realistic experiments, compared to the single-species case single_component .

One typically would not expect such supersolids as the bare species alone can only form superfluid and MI phases due to the short range two-body interactions two_component_SF_MI . The underlying mechanism is that the flow of the bare species is suppressed by the crystalline distribution of the dressed species via the interspecies interaction. As a result, the widths of the SS1 and SS2 phases will strongly depend on the interspecies interaction $U_{\text{bd}}$. The numerical result in Fig. 3b shows that indeed the two SS phases shrink as $U_{bd}$ decreases. The SS1 phase eventually disappears for sufficiently small $U_{\text{bd}}$.

For even larger hopping rate $t$, both species are in SF phases, which are characterized by nonzero SF order parameters. Different from the SS, spatial densities of both species become homogeneous in the SF phases.

Supersolidity mechanism of the bare species— In the rest of the work, we will develop a Bogoliubov mean-field theory to understand how the interspecies interaction enables the bare species to form SS phases. Our discussion will focus on the SS1 phase, where the dressed species is a DW. This allows us to write down wave functions $|{\rm{DW}}_{d}\rangle$ of the DW according to the crystalline structure. We also assume that the total wave function in the ground state can be decoupled as $|\Psi_{\rm{g}}\rangle\approx|{\rm{DW}}_{d}\rangle\otimes|\Psi_{b}\rangle$, where $|\Psi_{b}\rangle$ is the wave function of the bare component. Then we can derive an effective Hamiltonian for the bare species by tracing out the dressed atom part, i.e. $\hat{H}_{\rm{e}}=\langle{\rm{DW}}_{d}|\hat{H}|{\rm{DW}}_{d}\rangle$. Explicitly the effective Hamiltonian reads,

[TABLE]

where ${\{j\}}$ denotes lattice sites occupied by dressed atoms. For convenience, we have omitted the index $b$ of the bare species. The last term gives the interspecies interaction, where the mean particle number per site of the dressed atoms $n_{d}=1$ has been used explicitly. A constant term, $C=\langle{\rm{DW}}_{d}|\sum_{i<j}V_{ij}\hat{n}_{id}\hat{n}_{jd}|{\rm{DW}}_{d}\rangle$ characterizing the long-range interaction energy, is neglected in the effective Hamiltonian.

The interaction with the dressed atoms (the last term in the effective Hamiltonian) introduces a new spatially periodic structure to the bare species, in addition to the optical lattice. As an example, we consider parameters corresponding to Fig. 2b. Here, the dressed atoms form an oblique lattice, see Fig. 3c for a cartoon picture of the 2D structure. The primitive cell of the new oblique lattice is apparently larger than the original lattice. In this example, the primitive lattice vectors are $a_{1}=(1,\,4)$ and $a_{2}=(4,\,1)$, with which we obtain the area of the primitive lattice $A=|a_{1}\times a_{2}|=15$, while the area of the optical lattice is 1. In turn, the corresponding reciprocal lattice is smaller than that of the optical lattice. To illustrate this, we plot the first Brillouin zone of the two lattices in Fig. 3b. Apparently they overlap only in a small central area (low momentum regions).

As a result, phonon excitations for momentum components in and out of the overlap region will be very different. To show this, we calculate the Bogoliubov dispersion relation of the effective Hamiltonian. In the low momentum region (where the two Brillouin zones overlap), $E_{k}^{\text{i}}=\sqrt{\varepsilon_{k}^{2}+2\bar{n}_{b}U\varepsilon_{k}}$ with $\varepsilon_{k}=-2t(\cos k_{x}a+\cos k_{y}a-2)$. Outside this region, the dispersion becomes $E_{k}^{\text{o}}=\sqrt{\left(\varepsilon_{k}-\bar{n}_{d}U_{bd}\right)^{2}+2\bar{n}_{b}U\left(\varepsilon_{k}-\bar{n}_{d}U_{bd}\right)}$. Here $\bar{n}_{b}$ ($\bar{n}_{d}$) are the mean population of the bare (dressed) component. Consequently, the dispersion is not continuous any more at the boundary of the Brillouin zone of the oblique lattice. The dispersion relation becomes complex when $U_{bd}>2t\bar{n}_{d}(2-\cos k^{(b)}_{x}a-\cos k^{(b)}_{y}a)$ where $k^{(b)}_{x}$ and $k^{(b)}_{y}$ are momenta at the boundary. In Fig. 1d, we plot the dispersion relation along the $k_{x}$ axis by varying the interspecies interaction $U_{bd}$, where the mode frequency becomes complex at $U_{bd}=U$. This so-called roton-like instability SF_SS here indicates that the emergence of supersolids is indeed induced by the strong interspecies interaction. Note that the mechanism here is different from SS phases induced by geometrically dependent hopping found in frustrated lattices frustration .

Interaction potentials of Rydberg dressed atoms— The level structure used in the Rydberg dressing is shown in Fig. 1a. The species $|d\rangle$ is coupled to a Rydberg state by an off-resonant laser with Rabi frequency $\Omega$ and detuning $\Delta$. Interactions between Rydberg atoms are of van der Waals type $V_{\text{r}}=C_{6}/r^{6}$, where $C_{6}$ is the respective dispersion coefficient. The Rydberg dressing gives the soft-core interaction $V_{ij}$ where the effective dispersion coefficient $\tilde{C}_{6}=(\Omega/\Delta)^{4}C_{6}$ and soft-core radius $R_{c}=(C_{6}/2\Delta)^{1/6}$. $R_{c}$ varies with the Rydberg states and detuning. For example, one can choose the Rydberg 36S state of 87Rb atoms ($C_{6}=241.6\,\text{MHz}\times\mu m^{6}$) and lattice constant $a=532$ nm. When $\Delta=7$ MHz, we obtain $R_{c}\approx 3a$. With this fixed detuning $\Delta$, the strength of the soft-core interaction is now controlled by the Rabi frequency $\Omega$.

To probe different phases shown in Fig. 2, one needs to change the parameters $V$, $U$ and $t$ together or separately over certain ranges. One simple way to achieve this is to tune the lattice potential depth $V_{0}/E_{r}$. In optical lattices, the onsite interaction $U$ depends on the lattice depth through $U=\sqrt{8/\pi}ka_{s}E_{r}(V_{0}/E_{r})^{3/4}$ and the hopping rate $t$ through $t=4/\sqrt{\pi}E_{r}(V_{0}/E_{r})^{3/4}\exp[-2(V_{0}/E_{r})^{1/2}]$ interaction_hopping , where $k=2\pi/\lambda$, $E_{r}=h^{2}/2m\lambda^{2}$, $\lambda=2a$ and $a_{s}$ are the wave number, recoil energy, wavelength of the lattice potential and s-wave scattering length, respectively. Upon varying $V_{0}/E_{r}$ and fixing the other parameters, the ratios $t/U$ and $V/U$ change continuously. One example is shown in Fig. 4. One can see that the parameters cross the main phases discussed in this paper.

In conclusion, we have investigated crystalline phases of ultracold binary bosonic gases on a square lattice, with one species possessing a non-local interaction induced by Rydberg dressing. We found two types of supersolid phases that are robust and occupy large parameter regions at zero temperature. We showed that the supersolid phases of the bare species are stabilized by the interspecies interaction. The existence of the different phases predicted here could be directly observed by quantum gas microscopy with single-site resolution single_site ; M. Greiner_2009 ; I. Bloch_2010 or through measuring noise correlations noise . Our results demonstrate rich features of the Bose-Bose mixture with long-range interactions, and indicate that this system is well suited for exploring supersolidity in upcoming experiments. As the crystalline structure (see Fig. 3a) can be changed in the insulating region by tuning $V/U$, we expect that supersolid phases with tunable density patterns can be explored as well.

We acknowledge J.-M. Yuan, Z.-X. Zhao and Rejish Nath for useful discussions. This work was supported by the National Natural Science Foundation of China under Grants No. 11304386 and No. 11104350. WL acknowledges support from the UKIERI-UGC Thematic Partnership No. IND/CONT/G/16-17/73 and EPSRC Grant No. EP/M014266/1. AG and WH acknowledge support from the Deutsche Forschungsgemeinschaft via DFG SPP 1929 GiRyd and the high-performance computing center LOEWE-CSC.

I Supplementary Material

II Method

II.1 RBDMFT equations

In deriving the effective action, we consider the limit of a high but finite dimensional optical lattice, and use the cavity method georges96 ; Byczuk_2008 to derive self-consistency equations within RBDMFT. In a more formal language, first we map the Hamiltonian onto a set of individual single-site problems each of which is described by a local effective action LYQ :

[TABLE]

Here we have defined the local Weiss Green’s function,

[TABLE]

and introduced

[TABLE]

as the superfluid order parameters, and

[TABLE]

as the diagonal and off-diagonal parts of the connected Green’s functions, respectively, where $\langle\ldots\rangle_{0}$ denotes the expectation value in the cavity system (without the impurity site) LYQ ; Walter .

II.2 Anderson impurity model

The most difficult step in the procedure discussed above is to find a solver for the effective action. However, one cannot do this analytically. To obtain RBDMFT equations, it is better to return back to the Hamiltonian representation. Here, each of the local effective actions (II.1) is represented by an Anderson impurity Hamiltonian

[TABLE]

where the chemical potential and interaction term are directly inherited from the Hubbard Hamiltonian. The bath of condensed bosons is represented by the Gutzwiller term with superfluid order parameters $\phi_{\sigma}$ for each component. The bath of normal bosons is described by a finite number of orbitals with creation operators $\hat{a}^{\dagger}_{l}$ and energies $\epsilon_{l}$, where these orbitals are coupled to the impurity via normal-hopping amplitudes $V_{\sigma,l}$ and anomalous-hopping amplitudes $W_{\sigma,l}$. The anomalous hopping terms are needed to generate the off-diagonal elements of the hybridization function. Note here that in the high-dimensional limit inter-site interactions only contribute to the Hartree level Hartree . In other words, the Hartree term of the inter-site interaction will dominate as the spatial dimension of the system increases. This motivates us to keep only the Hartree contribution of the inter-site interaction in our simulations as an approximation to the original Hamiltonian, i.e.

[TABLE]

We now turn to the solution of the impurity model. In practice, we start with an initial set of Anderson paramters and local bosonic superfluid order parameters $\phi_{j,\nu}(\tau)$. The Anderson Hamiltonian can straightforwardly be implemented in the Fock basis, and the corresponding solution can be achieved by exact diagonalization (ED) of DMFT M. Caffarel_1994 ; georges96 . After diagonalization, the local Green’s function, which includes all the information about the bath, can be obtained from the eigenstates and eigenenergies in the Lehmann-representation

[TABLE]

Integrating out the orbitals leads to the same effective action as in Eq. (II.1), if the following identification is made

[TABLE]

where $\mathbf{G}_{\sigma\sigma^{\prime},ij}(i\omega_{n})$ is the inverse Fourier transformation of the Weiss Green’s function defined in Eq. (4) and (5), and the hybridization functions read:

[TABLE]

Hence, we obtain a set of local self-energies $\Sigma^{(i)}_{{\rm imp},\sigma\sigma^{\prime}}(i\omega_{n})$,

[TABLE]

Then we employ the Dyson equation in real-space representation in order to compute the interacting lattice Green’s function

[TABLE]

The site-dependence of the Green’s functions is shown by boldface quantities that denote a matrix form with site-indexed elements. Here $\mathbf{G}_{0}(i\omega_{n})^{-1}$ stands for the inverse non-interacting Green’s function

[TABLE]

In this expression, $\mathbf{1}$ is the unit matrix, the matrix elements $t_{ij}$ are hopping amplitudes for a given lattice structure. Eventually the self-consistency loop is closed by specifying the Weiss Green’s function via the local Dyson equation

[TABLE]

where the diagonal elements of the lattice Green’s function yield the interacting local Green’s function $\mathbf{G}^{(i)}_{\sigma\sigma^{\prime}}(i\omega_{n})=(\mathbf{G}_{\sigma,\sigma^{\prime}}(i\omega_{n}))_{ii}$. This self-consistency loop is repeated until the desired accuracy for superfluid order parameters and Anderson parameters is obtained.

II.3 Energy within RBDMFT

Calculation of energy is not straightforward within RBDMFT, since the kinetic energy kinetic is given in terms of non-local expectation values. It can be shown that within the RBDMFT self-consistency conditions, kinetic energy can also be written in terms of Anderson impurity hybridization functions and local Green’s functions. A detailed derivation can be found in Ref. Huber .

II.3.1 Kinetic energy

In terms of creation and annihilation operators for bosons, $b^{\dagger}_{i\sigma}$ and $b_{i\sigma}$, respectively, kinetic energy has the form

[TABLE]

Thus expressing the total kinetic energy in terms of real-space Green’s functions yields

[TABLE]

This expression can be further simplified by employing both the local and lattice Dyson equations within RBDMFT

[TABLE]

which yields

[TABLE]

Further using the self-consistency property of the impurity Green’s function leads to

[TABLE]

and we finally obtain

[TABLE]

II.3.2 Total energy

The ground state within RBDMFT corresponds to the solution with the lowest energy, where the corresponding total energy of the impurity site which is given as follows:

[TABLE]

For the Bose-Hubbard model of spin-1 bosons, the on-site interaction term is given by:

[TABLE]

III Numerical results within RBDMFT

III.1 Density dependent phase diagram of Rydberg-dressed systems

In this paragraph, we study the stability of quantum phases of Rydbery-dressed systems in optical lattices for different fillings. In the strong coupling limit with $U_{\sigma\sigma^{\prime}}\gg t$, we find that the system favors Mott insulating or density-wave phase with different types of crystalline order in the individual species. Interestingly, we observe a density-wave phase with a nonuniform total density which breaks lattice translational symmetry, with densities $n_{ib}=1$ and $n_{id}=2$, appears, as shown in green region of Fig. S1. These density waves exhibit nonzero density fluctuations, as shown in Fig. S2. However, quantum fluctuations as a result of higher-order tunneling processes are weak, due to the strong long-range interactions. Actually, the density wave of the dressed species is also predicted in the single-species case long-range .

Away from the deep MI regime, i.e. in the intermediate hopping regime, we observe two types of quantum phase transition from MI to supersolid, i.e. the uncoupled ground-state species demonstrates a phase transition from MI to supersolid, and then followed by the Rydberg dressed species, as shown in the Fig. S1. Interestingly, we observe a pronounced region of supersolid appearing in our simulations, as a result of the onsite interspecies interactions, indicating a higher chance for directly observing these phases in realistic experiments, compared to single-species case long-range . Actually, we indeed observe the width of SS1 and SS2 shrinks as a function of interspecies interactions, as shown in Fig. S1(f), where SS1 clearly disappears for smaller $U_{bd}$. In addition, the long-range interaction also shifts the phase transition between MI and SS1, even though the bare species only possess onsite interactions. As shown in Fig. 2 in the main text, the phase boundary shrinks to lower hopping regime with increasing the long-range interaction $V$.

Finally, in the weakly interacting regime with $t\gg U_{\sigma\sigma^{\prime}}$, a superfluid phase with uniform total density distribution is found in our simulations, where both species demonstrate homogenous density distribution. Here, crystalline orders are destroyed by the large density fluctuations, and the system only supports superfluidity with uniform density.

III.2 Dipolar system

We have so far studied crystalline order in the Rydberg dressed systems. Actually, the physics of these competing orders can also be exhibited in dipolar system loaded in an optical lattice, along with quick developments in the cooling and trapping of magnetic atoms Cy_Dr_Ey and diatomic molecules Ye . Recently, a Gutzwiller mean-field phase diagram of a binary Bose mixture on a square optical lattice is studied, where one species possesses a non-negligible dipole moment NN . In their study, only the nearest-neighbor part of the dipolar interactions was included. To obtain a better understanding of the Rydberg dressed system studied above and make a comparison, we here study a mixture of dipolar and nondipolar bosons on a square optical lattice, with real long-range interactions beyond nearest-neighbor approximations. We study the system by means of RBDMFT, which takes into account quantum fluctuations and is actually a higher-order expansions of Gutzwiller mean-field theory.

In Fig. S3, we show the resulting phase diagram of dipolar and nondipolar bosonic mixtures on a 2D optical lattice. In general, there are also five phases in this dipolar system, i.e. SF, MI, DW and two types of supersolid. Compared to nearest-neighbor interaction and static mean-field approximations NN , two big differences have been observed. First, rich crystalline patterns appear in the system, as shown in Fig. S3a)-c), with a filling factor of 1/3, 1/4 and 1/8 for the dipolar species, respectively. Second, we observe that the region of supersolid phase is also altered. Note here that we recover the static mean-field phase diagram with nearest-neighbor interactions within Gutzwiller approximations in Ref. NN .

IV Bogoliubov spectra of the bare species in the SS1 phase

As the dressed atoms are in a density wave state, we could decouple the total wave function in the ground state as $|\Psi_{\rm{g}}\rangle\approx|{\rm{DW}}_{d}\rangle\otimes|\Psi_{b}\rangle$, where $|{\rm{DW}}_{d}\rangle$ and $|\Psi_{b}\rangle$ are the wave function of dressed atoms and bare component, respectively. Here quantum fluctuations of the density wave could be neglected. After tracing out the dressed atoms, we obtain an effective Hamiltonian for the bare species,

[TABLE]

where ${\{j\}}$ denotes sites of the oblique lattice occupied by dressed atoms with the corresponding particle number $n_{d}$. For parameters considered in this work, numerical results show that $n_{d}\approx 1$. Through Fourier transformation, we can derive the Hamiltonian in momentum space in the first Brillouin zone,

[TABLE]

where $N$ is the total number of sites and $U_{1}=\bar{n}_{d}U_{bd}$ with $\bar{n}_{d}=n_{d}N_{d}/N$. $N_{d}$ is the number of sites occupied by the dressed atoms, and $\{k\}$ denotes momentum spanned in the first Brillouin zone of the lattice occupied by the dressed atoms.

Expanding the Hamiltonian (S32) around $|\vec{k}|=0$ and keeping only quadratic terms of the operators, this yields,

[TABLE]

where $E_{0}=-UN_{0}^{2}/2N$ is the energy of the condensed atoms, with $N_{0}$ to be the number of condensed atoms and $\mu=-4t+U\bar{n}_{b}+U_{1}$ the chemical potential and the mean occupation of the condensed atom $\bar{n}_{b}=N_{0}/N$.

As the interspecies interaction [the last term in Eq. (S33)] only appears in the low momentum regions (Brillouin zone $\{k\}$), we will have two different forms of the approximate Hamiltonian depending on values of the momentum. Substituting the chemical potential $\mu$, we get the approximate Hamiltonian within the first Brillouin zone of the dressed atom,

[TABLE]

and the corresponding Bogoliubov spectrum is

[TABLE]

with $\varepsilon_{k}=-2t(\cos k_{x}a+\cos k_{y}a-2)$. The spectrum is similar to the one of a weakly interacting Bose gas in a square optical lattices.

For momenta outside the first Brillouin zone of the dressed atoms, we have a different form of the approximate Hamiltonian,

[TABLE]

the corresponding Bogoliubov spectra is

[TABLE]

which will be nonzero only at large momentum (outside the first Brillouin zone).

The roton instability occurs at the boundary of the two Bogoliubov spectrum. Using Eq. (S37), we can find the spectrum becomes complex when $\varepsilon_{k}<U_{1}$. This allows us to find the critical value of the tunneling rate $t_{c}$

[TABLE]

where $k^{(b)}_{x}$ and $k_{y}^{(b)}$ are values of the momentum at the boundary of the first Brillouin zone of the oblique lattice.

The soft-core interaction will affect structures of the oblique lattice. Therefore the critical $t_{c}$ will change as the interaction $V$ changes. As shown in Fig. 3 in the main text, the first Brillouin zone is not of a regular shape, such that the critical value $t_{c}$ will vary with both $k^{(b)}_{x}$ and $k^{(b)}_{y}$. To show this we evaluate the critical values $t_{c}$ using the crystalline structure of the dressed atoms at the SS1-SS2 phase boundary, which are obtained by the full numerical calculation. For example, $t_{c}$ lies in a range $[0.087,0.094]$ if $V=0.3$. When the long range interaction becomes strong, we find that the range of critical $t_{c}$ increases. For example, $t_{c}\in[0.085,0.11]$ when $V=0.4$, and $t_{c}\in[0.073,0.13]$ when $V=0.6$. Although these values are close to the numerical calculations, it is apparent that one will not be able to determine phase boundaries accurately using the Bogoliubov calculation.

Another limitation of this calculation is that areas of the crystalline structure become smaller when $V$ is weak. Long range correlations become important in determining the ground state phases, which prevents us to decouple the total wave function into two parts. In this regime, the Bogoliubov calculation fails to capture the many body physics.