Mean-field phase diagram of ultracold atomic gases in cavity quantum electrodynamics

TL;DR

This paper maps the phase diagram of ultracold bosonic atoms in a 2D optical lattice coupled to a cavity, revealing the stability of supersolid phases and comparing with previous theoretical and numerical studies.

Contribution

It provides a comprehensive mean-field analysis of the phase diagram for a 2D Bose-Hubbard model with infinite-range interactions, including analytical and numerical results.

Findings

Agreement with previous phase diagrams from literature

Identification of stability regions for supersolid phases

Discussion of differences with dipolar boson systems

Abstract

We investigate the mean-field phase diagram of the Bose-Hubbard model with infinite-range interactions in two dimensions. This model describes ultracold bosonic atoms confined by a two-dimensional optical lattice and dispersively coupled to a cavity mode with the same wavelength as the lattice. We determine the ground-state phase diagram for a grand-canonical ensemble by means of analytical and numerical methods. Our results mostly agree with the ones reported in Dogra et al. [PRA 94, 023632 (2016)], and have a remarkable qualitative agreement with the quantum Monte Carlo phase diagrams of Flottat et al. [PRB 95, 144501 (2017)]. The salient differences concern the stability of the supersolid phases, which we discuss in detail. Finally, we discuss differences and analogies between the ground state properties of strong long-range interacting bosons with the ones predicted for repulsively interacting dipolar bosons in two dimensions.

Figures

Equations

$${\hat{H}}_t$$

$${\hat{V}}_0$$

$$\hat{\Phi }=\sum _j(-1{)}^j{\hat{n}}_j/K,$$

$${\hat{H}}_{GC}=\hat{H}-\mu \sum _j{\hat{n}}_j.$$

$$\rho =\frac{1}{K}\sum _j⟨{\hat{n}}_j⟩,$$

$$\theta =2\left|⟨\hat{\Phi }⟩\right|,$$

$$n_j=\rho +(-1{)}^j\frac{\theta }{2}+{\delta }_j,$$

$$\sum _j{\delta }_j=0,$$

$$\sum _j(-1{)}^j{\delta }_j=0.$$

$$E(\rho ,\theta ,\{{\delta }_j\})=E_0(\rho ,\theta )+\frac{U_0}{2}\sum _j{\delta }_j^2,$$

$$E_0(\rho ,\theta )=K\left[\frac{U_0}{2}\rho (\rho -1)+\left(\frac{U_0}{2}-U_{\infty }\right)\frac{{\theta }^2}{4}-\mu \rho \right],$$

$$U_0(n-1)+\frac{U_{\infty }}{2}<\mu <U_{0}n-\frac{U_{\infty }}{2}.$$

$$(U_0-U_{\infty })(\theta -1)-\frac{U_{\infty }}{2}<\mu <(U_0-U_{\infty })\theta -\frac{U_{\infty }}{2},$$

$${\phi }_i=⟨{\hat{a}}_i⟩.$$

$${\hat{H}}_t^{(i)}=-t\left({\hat{a}}_i^{†}{\bar{\phi }}_i+{\hat{a}}_i{\bar{\phi }}_i^{\ast }-\mathrm{Re}\{{\phi }_i^{\ast }{\bar{\phi }}_i\}\right),$$

$${\hat{H}}_{\mathrm{mf}}^{(i)}={\hat{H}}_t^{(i)}+\frac{U_0}{2}{n}_i(n_i-1)-(-1{)}^i{U}_{\infty }\theta n_i+U_{\infty }\frac{{\theta }^2}{4}-\mu n_i.$$

$$\begin{array}{rl}{H}_{s\in \{\text{e},\text{o}\}}=& -zt{\phi }_{\bar{s}}({\hat{a}}_s+{\hat{a}}_s^{†}-{\phi }_s)+\frac{U_0}{2}{\hat{n}}_s({\hat{n}}_s-1)\\ & -\mu {\hat{n}}_s-{\sigma }_s{U}_{\infty }\theta {\hat{n}}_s+U_{\infty }{\theta }^2/4\text{,}\end{array}$$

$$Z=\mathrm{Tr}\left\{e^{-\beta {\hat{H}}_0}{T}_{\tau }{e}^{-{\int }_0^{\beta }d\tau {\hat{H}}_I(\tau )}\right\}$$

$${\hat{H}}_I(\tau )=e^{\tau {\hat{H}}_0}{\hat{H}}_t{e}^{-\tau {\hat{H}}_0},$$

$$Z=\int \mathcal{D}{\alpha }_j\mathcal{D}{\alpha }_j^{\ast }{e}^{-{\int }_0^{\beta }d\tau \mathcal{L}(\tau )}$$

$$\mathcal{L}=\sum _j{\alpha }_j^{\ast }\frac{d{\alpha }_j}{d\tau }+H(\{{\alpha }_j^{\ast },{\alpha }_j\}).$$

$$\begin{gathered} ⟨T_{\tau }A[\{{\hat{a}}_j^{†}({\tau }_j),{\hat{a}}_{j^{\prime }}({\tau }_{j^{\prime }})\}]{⟩}_0= \\ \frac{1}{Z_0}\int \mathcal{D}{\alpha }_j\mathcal{D}{\alpha }_j^{\ast }{e}^{-{\int }_0^{\beta }d\tau \mathcal{L}(\tau )}A[\{{\alpha }_{j^{\prime }}^{\ast }({\tau }_{j^{\prime }}),{\alpha }_j({\tau }_j)\}], \end{gathered}$$

$${\alpha }_q=\frac{1}{\sqrt{K}}\sum _j{\alpha }_j\exp [2\pi i(j_1{q}_1+j_2{q}_2)/\sqrt{K}]$$

$$\mathcal{L}={\mathcal{L}}_0-t\sum _q{w}_q{\alpha }_q^{\ast }{\alpha }_q,$$

$$e^{t\int d\tau {\sum }_q{w}_q{\alpha }_q^{\ast }{\alpha }_q}=\int \mathcal{D}{\psi }_q\mathcal{D}{\psi }_q^{\ast }{e}^{-\int d\tau {\mathcal{L}}_2+{\mathcal{L}}_c},$$

$${\mathcal{L}}_2$$

$${\mathcal{L}}_c$$

$$Z=Z_0\int \mathcal{D}{\psi }_q\mathcal{D}{\psi }_q^{\ast }{e}^{-S_{\mathrm{eff}}},$$

$$S_{\mathrm{eff}}$$

$$S_{\rm eff}^{(2)}=-\frac{1}{2}\Big{\langle}\Big{(}\int_{0}^{\beta}d\tau\mathcal{L}_{c}\Big{)}^{2}\Big{\rangle}_{0}+\int_{0}^{\beta}d\tau\mathcal{L}_{2}$$

$$S_{\rm eff}^{(2)}=\int_{0}^{\beta}\!\!\!d\tau\mathcal{L}_{2}-t^{2}\sum_{q,q^{\prime}}\sqrt{w_{q}w_{q^{\prime}}}\int_{0}^{\beta}\!\!\!d\tau\int_{0}^{\tau}\!\!\!d\tau^{\prime}\\ \Big{[}\psi_{q}^{*}(\tau)\psi_{q^{\prime}}(\tau^{\prime})\langle\hat{a}_{q}(\tau)\hat{a}_{q^{\prime}}^{\dagger}(\tau^{\prime})\rangle_{0}+\psi_{q}(\tau)\psi_{q^{\prime}}^{*}(\tau^{\prime})\langle\hat{a}^{\dagger}_{q}(\tau)\hat{a}_{q^{\prime}}(\tau^{\prime})\rangle_{0}\Big{]}\,,$$

$$⟨{\hat{a}}_j^{†}(\tau ){\hat{a}}_{j^{\prime }}(\tau -{\tau }_0){⟩}_0$$

Full text

I Introduction

The Bose-Hubbard model is a paradigmatic quantum mechanical description of strongly-correlated spinless particles in a lattice bloch-dalibard-zwerger-2008 . This model predicts a quantum phase transition, which emerges from the competition between the hopping coupling nearest-neighbour lattice sites and the onsite repulsion. The latter penalizes multiple occupation at the individual sites and favours the Mott-insulator (MI) phase, while the hopping promotes the buildup of non-local correlations and the onset of superfluidity (SF) fisher-et-al-1989 . Ultracold vapours of alkali-metal atoms offer a remarkable setup for realising and investigating the dynamics of the Bose-Hubbard model, thanks to the possibility to independently tune onsite interactions and lattice depth jaksch-et-al-1998 ; greiner-et-al-2002 . The addition of further interactions, such as interparticle potentials that decay with a power-law of the inverse of the distance, typically gives rise to frustration and results in the appearance of new phases menotti-et-al-2008 ; li-et-al-2018 ; larson-et-al-2008 . The recent realization and measurement of the phases of ultracold atomic lattices in an optical cavity has set the stage for the studies of Bose-Hubbard models with infinite-range interactions landig-et-al-2016 . Here, attractive infinite-range interactions are realised via multiple scattering of cavity photons and favour the formation of density patterns which maximize the intracavity field domokos-ritsch-2002 ; asboth-et-al-2005 ; fernandez-vidal-et-al-2010 ; habibian-et-al-2013 . When the atoms are confined by an external optical lattice and the wavelengths of the cavity mode and of the laser generating the underlying optical lattice coincide, the emerging patterns can have a checkerboard density modulation, and are denoted by charge-density wave (CDW) or lattice supersolid (SS) depending on whether the gas is incompressible or superfluid, respectively.

The experimentally measured phase diagram of Refs. landig-et-al-2016 ; hruby-et-al-2018 reports the existence of these four phases. Several theoretical works reproduced the salient features of the phase diagram using different approaches. Most works use different implementations of the mean-field treatment dogra-et-al-2016 ; niederle-morigi-rieger-2016 ; sundar-mueller-2016 ; wald-et-al-2018 ; liao-et-al-2018 , nevertheless their predictions do not agree across the whole phase diagram. Moreover, they exhibit several discrepancies with state-of-the-art two-dimensional quantum Monte Carlo (QMC) study flottat-et-al-2017 . It has been further argued that the mean field predictions for this kind of Bose-Hubbard model shall be the same as the one of Bose-Hubbard with nearest-neighbour (and thus also repulsive dipolar) coupling dogra-et-al-2016 ; niederle-morigi-rieger-2016 .

The objective of this work is to provide an extensive mean-field analysis of the two-dimensional Bose-Hubbard model with infinite long-range interactions, describing the setup of Ref. landig-et-al-2016 . For this analysis we consider a grand-canonical Hamiltonian and determine the ground-state phase diagram by taking particular care of the convergence criterion of the numerical results and by comparing them with analytical results. We then discuss our results comparing them in detail with the analytical predictions and with previous theoretical studies for the cavity Bose-Hubbard model and for repulsively interacting dipolar gases in two dimensions.

This manuscript is organized as follows. In Sec. II we introduce the Bose-Hubbard model and review the ground-state properties in the atomic limit, namely, when the kinetic energy is set to zero. In Section III we derive the mean-field Hamiltonian and employ the path-integral formalism to analytically determine the transition from incompressible to compressible phases. In Sec. IV we numerically determine the ground-state phase diagram and compare our results with the ones reported so far in the literature. The conclusions are drawn in Sec. V while the appendices provide details on the numerical methods for calculating the mean-field phase diagrams.

II Bose-Hubbard Hamiltonian with cavity-mediated interactions

In this section we introduce the Bose-Hubbard model in presence of cavity-mediated long-range interactions and review the exact result for eigenstates and eigenspectrum in the so-called atomic limit, where the kinetic energy is set to zero. This limit is relevant for the mean-field study, where we determine the ground state for finite values of the hopping term.

II.1 Grand-canonical Hamiltonian

We consider $N$ bosons tightly confined in the lowest band of a two-dimensional, quadratic optical lattice. The lattice has $K=L\times L$ sites, each site is labeled by the index $i=(i_{1},i_{2})$ and we assume periodic boundary conditions. The bosons are described by operators $\hat{a}_{i}$ and $\hat{a}_{i}^{\dagger}$, which annihilate and create a particle at site $i=(i_{1},i_{2})$, respectively, and obey the commutation relation $[\hat{a}_{i},\hat{a}_{j}^{\dagger}]=\delta_{i_{1},j_{1}}\delta_{i_{2},j_{2}}$. The bosons interact via $s$-wave scattering. They also interact via the long-range interactions mediated by the photons of a single-mode cavity, whose periodicity is commensurate with the optical lattice and which gives rise to all-to-all coupling. The Hamiltonian $\hat{H}=\hat{H}_{t}+\hat{V}_{0}$ governing the many-body dynamics is here decomposed into the kinetic energy $H_{t}$ and potential energy $V_{0}$, which individually read habibian-et-al-2013

[TABLE]

The kinetic energy is scaled by the hopping coefficient $t$, which is positive and uniform across the lattice, and $\sum_{\langle ij\rangle}$ is restricted to the pairs of nearest neighbour sites $i$ and $j$. The potential energy is diagonal on the eigenstates of operator $\hat{n}_{i}=\hat{a}_{i}^{\dagger}\hat{a}_{i}$ and consists of the onsite repulsion, which is scaled by the strength $U_{0}>0$, and of the attractive infinite-range interactions with strength $U_{\infty}$, which multiplies the extensive operator $K\hat{\Phi}^{2}$. For the setup of Ref. landig-et-al-2016 , where cavity and optical-lattice laser wave lengths are equal, operator $\hat{\Phi}$ takes the form

[TABLE]

where $(-1)^{j}\equiv(-1)^{j_{1}+j_{2}}$. Its expectation value is maximum and proportional to the mean density when the atoms form a checkerboard pattern. We denote a site $(j_{1},j_{2})$ by even (odd) when $j_{1}+j_{2}$ is an even (odd) number.

In the rest of this paper we will study the phase diagram of a grand-canonical ensemble at $T=0$. For this purpose we analyse the ground state of the grand-canonical Hamiltonian, defined as

[TABLE]

Here, $\mu$ is the chemical potential which controls the mean occupation number $\rho$,

[TABLE]

and the expectation value is taken over the grand-canonical ensemble. In the following we also use the parameter $\theta$, which is proportional to the expectation value of operator $\hat{\Phi}$ according to the relation:

[TABLE]

where the proportionality factor 2 is introduced for later convenience. The value of $\theta$ measures the population imbalance between even odd sites, thus when it is non-vanishing the atomic density is spatially modulated. In particular, it is proportional to the value of the structure form factor at the wave number of the cavity field niederle-morigi-rieger-2016 .

II.2 Atomic limit

We now consider the limit $t=0$. In this case the energy eigenstates are the Fock states $|n_{(1,1)},\ldots,n_{(L,L)}\rangle$, with $|n_{i}\rangle$ Fock state at site $i$. It is convenient to decompose the Fock number $n_{j}$ of each site $j$ as the sum

[TABLE]

where $\delta_{j}$ ensures that $n_{j}$ is a natural number. This condition, together with Eqs. (5) and (6), lead to the relations

[TABLE]

Using these relations and Eq. (7) one can verify that, whenever $\rho\pm\theta/2$ is an integer number, the configuration with minimal energy has $\delta_{j}=0$. In fact, using Eq. (7) one can cast the energy of the state $|n_{(1,1)},\ldots,n_{(L,L)}\rangle$ into the form

[TABLE]

where $E_{0}(\rho,\theta)$ is the energy of the configuration when $\delta_{j}$ vanishes at all sites,

[TABLE]

and is visibly extensive. This expression is correct when $|\theta|\leq 2\rho$. The ground state is found by the configuration which minimizes the energy $E_{0}(\rho,\theta)$. Therefore the ground-state properties are determined by two independent parameters, which we choose here to be $U_{\infty}$ and $\mu$ in units of $U_{0}$. More generally, the states at energy $E_{0}$ are characterised by two-site translational symmetry along both directions of the lattice, such that the sites with the same parity are equally populated. Hence, we can denote the ground state by the ket $\{n,m\}$ where $n$ ($m$) is the Fock number for the even (odd) sites, or vice versa.

In the following we review the ground-state properties in the thermodynamic limit $K\to\infty$ by analysing Eq. (11). They can be displayed by means of a phase diagram in the $U_{\infty}-\mu$ space shown in Fig. 1, see also Ref. dogra-et-al-2016 ; liao-et-al-2018 ; sundar-mueller-2016 . We first notice that in the limit $U_{\infty}=0$ the phase is MI with commensurate density $\rho=n$ in the interval $\mu\in[U_{0}(n-1),U_{0}n]$, while for $\mu<0$ the density is $\rho=0$. At $\mu=U_{0}n$ there is an infinite degeneracy of SF phases with density continuously varying from $n$ to $n+1$. For increasing value of $U_{\infty}$, but $U_{\infty}<U_{0}/2$, the MI phase with commensurate density $n$ is the ground state for values of the chemical potential such that

[TABLE]

At the upper (lower) boundary there is an abrupt jump from the MI to a CDW phase with fractional density $n+1/2$ ($n-1/2$) and population imbalance $|\theta|=1$. In this CDW the occupation of two adjacent sites is $\{n,n+1\}$ ($\{n,n-1\}$), or vice versa, the CDW ground state being doubly degenerate. These boundaries are the lines depicted in Fig. 1. At $U_{\infty}=U_{0}/2$ there is a discontinuity: For $U_{\infty}>U_{0}/2$ the ground state at density $\rho$ is a CDW with the largest population imbalance $|\theta|=2\rho$ (where $2\rho$ is an integer) in the interval

[TABLE]

while at $\mu=2(U_{0}-U_{\infty})\rho-\frac{U_{\infty}}{2}$ there is an infinite manifold of SF states with density varying from $\rho$ to $\rho+1/2$. The corresponding phases and boundaries are shown in Fig. 1, they all converge to the same point at $U_{\infty}=U_{0}$. For $U_{\infty}\geq U_{0}$ the onsite energy is attractive, the energy is not bound from below and the grand-canonical ensemble becomes unstable.

This analysis reproduces the results of Ref. dogra-et-al-2016 ; liao-et-al-2018 ; sundar-mueller-2016 for $t=0$. In what follows we will study the phase diagram for finite tunneling rates by means of the mean-field analysis.

III Mean-field analysis

In this Section we review the mean-field model which is at the basis of our numerical calculations and the definition of the observables that identify the phases. We then use the path-integral formalism to analytically determine the transition from incompressible to SF phase.

III.1 Mean-field treatment

We first introduce a so-called ”local superfluid order parameter”, which is the expectation value of the annihilation operator $\hat{a}_{i}$ at site $i$:

[TABLE]

The mean-field approximation consists of neglecting terms in second-order in the fluctuations $\delta\hat{a}_{i}$ of the annihilation operator about $\varphi_{i}$, with $\delta\hat{a}_{i}=\hat{a}_{i}-\varphi_{i}$. With this approximation the Hamiltonian term (1) can be cast into the sum of local Hamiltonians $\hat{H}_{t}^{\text{mf}}=\sum_{i}\hat{H}_{t}^{(i)}$ with

[TABLE]

and where $\bar{\varphi}_{i}=\sum_{\langle j\rangle_{i}}\varphi_{j}$ is the sum of local SF order parameters of the neighbours of site $i$. Without loss of generality in the numerical calculations we assume that these parameters are real.

In order to write the total Hamiltonian in terms of local operators, we perform a second approximation by writing the cavity potential in the mean-field form: $\hat{\Phi}^{2}\approx\theta\hat{\Phi}-\theta^{2}/4$, thus we discard fluctuations of $\hat{\Phi}$ to second order. With this approximation we can now write the grand-canonical Hamiltonian, Eq. (4), in its mean field form $H_{\rm GC:mf}=\sum_{i}\hat{H}^{(i)}_{\rm mf}$, namely, as the sum of local-site Hamiltonians $\hat{H}^{(i)}_{\rm mf}$ that read

[TABLE]

In the following we assume two-site symmetry, as in Ref. dogra-et-al-2016 . Using this assumption all even and all odd sites possess the energy $\hat{H}^{(e)}_{\rm mf}$ and $\hat{H}^{(o)}_{\rm mf}$, respectively, such that $\hat{H}_{\rm GC:mf}=K(\hat{H}^{(e)}_{\rm mf}+\hat{H}^{(o)}_{\rm mf})/2$. It is now convenient to introduce the annihilation and creation operators $\hat{a}_{e}$ and $\hat{a}_{e}^{\dagger}$ ($\hat{a}_{o}$ and $\hat{a}_{o}^{\dagger}$) for a particle in an even (odd) site, and the corresponding number operator $\hat{n}_{e}=\hat{a}_{e}^{\dagger}\hat{a}_{e}$ ($\hat{n}_{o}=\hat{a}_{o}^{\dagger}\hat{a}_{o}$). The even (odd) sites have SF order parameter $\varphi_{e}$ ($\varphi_{o}$) and the population imbalance operator reads $\hat{\Phi}=(\hat{n}_{e}-\hat{n}_{o})/2$. With these definitions we write

[TABLE]

where we have used that $\bar{\varphi}_{s}=z\varphi_{\bar{s}}$, with $z$ the coordination number (here equal to 4) and $\varphi_{\bar{e}}=\varphi_{o}$ ($\varphi_{\bar{o}}=\varphi_{e}$). Moreover, we have introduced the symbol $\sigma_{\text{e}}=+1$, $\sigma_{\text{o}}=-1$. Hamiltonian (17) is at the basis of the numerical results presented in the next section.

III.2 Transition from incompressible to compressible phases

We now determine the critical tunneling rate which separates compressible from incompressible phases. For this purpose we start from Hamiltonian (4) and consider an elaborate form of mean-field treatment following Refs. kampf-zimanyi-1993 ; fisher-et-al-1989 ; sachdev-2011 . We consider the partition function muehlschlegel-1978 ; fisher-et-al-1989 :

[TABLE]

where $\beta$ is the inverse temperature, $\tau$ is the imaginary time, $T_{\tau}$ is the imaginary-time ordering operator, $\hat{H}_{0}=\hat{V}_{0}-\mu\sum_{i}\hat{n}_{i}$ is the grand-canonical Hamiltonian without the kinetic energy, and we take $\hbar=1$ to simplify the notation. Moreover,

[TABLE]

where $\hat{H}_{t}$ is the tunneling Hamiltonian, Eq. (1). We can also write Equation (18) as $Z=Z_{0}\left\langle T_{\tau}e^{-\int_{0}^{\beta}d\tau\hat{H}_{1}(\tau)}\right\rangle_{0}$ fisher-et-al-1989 , where $Z_{0}$ the partition function for the model corresponding to $\hat{H}_{0}$ and the expectation value evaluated for the thermal state of the same model at inverse temperature $\beta$. Equivalently, one can cast the expression in terms of coherent-state path integrals kampf-zimanyi-1993 ; sachdev-2011 ; negele-orland-1998 :

[TABLE]

where

[TABLE]

where $H$ is assumed to be written in normal form and the path integral is over variables satisfying periodic boundary conditions. The two formalisms can be related by noting that for an operator $A[\{\hat{a}_{j}^{\dagger}(\tau_{j}),\hat{a}_{j^{\prime}}(\tau_{j^{\prime}})\}]$ negele-orland-1998 :

[TABLE]

where the imaginary time dependence of the operators is defined in the same way as in Eq. (19).

We define a new basis of Fourier-transformed variables $\alpha_{q},\alpha_{q}^{*}$, with $q=(q_{1},q_{2})$:

[TABLE]

with $L=\sqrt{K}$. We then write

[TABLE]

where $w_{q}=2(\cos{(2\pi q_{1}/\sqrt{K})}+\cos(2\pi q_{2}/\sqrt{K}))$ are the eigenvalues of the vicinity matrix, and where $\mathcal{L}_{0}$ is the Lagrangian without the tunneling terms. By means of the Hubbard-Stratonovich transformation we can write

[TABLE]

where all normalization factors are now included in the definition of the functional integral, and

[TABLE]

The prefactors here are chosen so that $\psi_{q}$ are dimensionless. In particular, the auxiliary variables $\psi_{q},\psi_{q}^{*}$ are related to the Fourier transform of the expectation values $\varphi_{i}$ by the equation $\langle\psi_{q}\rangle=\sqrt{w_{q}}\varphi_{q}$.

We now integrate Eq. (20) over the variables $\alpha_{j},\alpha_{j}^{*}$ and obtain

[TABLE]

where we have introduced the effective action $S_{\rm eff}$. The effective action is non-local in time and is given by the expression

[TABLE]

In order to find the transition points, it is sufficient to consider $S_{\rm eff}$ up to second order in the auxiliary fields. One recovers the expression kampf-zimanyi-1993 :

[TABLE]

Owing to the phase invariance of the model, Eq. (30) reduces to the form:

[TABLE]

where the Fourier-transformed operators $\hat{a}_{q}$ of the site operators $\hat{a}_{j}$ are defined in analogous form as Eq. (23).

The time correlators of the model with no hopping can be calculated easily in the site basis. For the case $T\to 0$ (i.e. $\beta\to\infty$) they are found to be:

[TABLE]

where the values of $n_{j}$ and $\langle\Phi\rangle$ are the ones that correspond to the ground state for $t=0$, see Sec. II. The energy $E_{j}^{\pm}$ is the energy variation resulting from the addition or subtraction of a particle at site $j$,

[TABLE]

where we neglected a term of order $1/K$ (note that $E_{j}^{-}$ is defined for $n_{j}>0$).

In the ground state, $n_{j}$ and $E_{j}^{\pm}$ only depend on the parity of the site, so one can cast the correlators of Eqs. (32) and (33) in the form:

[TABLE]

where the subindices correspond to $j$ being even or odd. This can be used to calculate the Fourier-transformed correlators:

[TABLE]

Here, we introduced the notation:

[TABLE]

and the sum of quasimomenta is taken to be mod $\sqrt{K}$. Thus, the presence of an even-odd asymmetry leads to non-vanishing correlators between momenta corresponding to $q,\bar{q}$. Note that for temperatures $T>0$ this structure is maintained, only the form of the single-site correlators is changed.

The correlators in Fourier basis can then be replaced in expression (31), and the sum can be made more compact by noting that $w_{\bar{q}}=-w_{q}$. We now make the standard assumption that the transition can be found by considering time-independent auxiliary fields $\psi_{q},\psi_{q}^{*}$, so that they can be taken out of the integrals. We obtain:

[TABLE]

For the case $T\simeq 0$, after performing the time integrals one gets:

[TABLE]

Thus, for each pair of modes $q,\bar{q}$, the effective action to second order has eigenvalues corresponding to a matrix of the form:

[TABLE]

with $\ell_{1},\ell_{2}$ the ($q$-independent) coefficients in Eq. (42). The smallest of each pair of eigenvalues reads $1-t|w_{q}|\sqrt{\ell_{1}^{2}-\ell_{2}^{2}}/2$. Hence, since the largest value of $|w_{q}|$ is equal to 4, the smallest eigenvalue of all pairs of modes in 2 dimensions is then $1-2t\sqrt{\ell_{1}^{2}-\ell_{2}^{2}}$. The transition point is found when this eigenvalue vanishes. After replacing the coefficients $\ell_{1},\ell_{2}$ one finds:

[TABLE]

This result coincides with the one reported in Refs. dogra-et-al-2016 ; sundar-mueller-2016 ; liao-et-al-2018 .

We conclude this section by remarking that this formalism should also allow one to identify the transition between SF and lattice SS. We have applied it in fact by treating the cavity potential as perturbation of the SF ground state and further performed a higher-order expansion of the effective action in order to analyse the effect of coupling between order parameters. The phase boundary we obtain, however, does not agree with the numerical results of the following section. For this reason we refrain to report the details of the calculation.

IV Ground-state phase diagram

In this section we determine the ground-state phase diagram using the mean-field model, Eq. (17). By rescaling the energy with $U_{0}$, the ground state is fully characterized by three parameters: $\mu$, which controls the density, $U_{\infty}$, which scales the cavity interactions, and the tunneling $t$. The numerical analysis is performed by identifying self-consistently the ground state using a fixed-point iteration detailed in Appendix A. By these means we identify four possible phases: (i) SF when $\varphi_{s}\neq 0$ and $\theta=0$; (ii) lattice supersolid (SS) when $\varphi_{s}\neq 0$ and $\theta\neq 0$; (iii) MI when $\varphi_{s}=0$ and $\theta=0$; and finally (iv) CDW, when $\varphi_{s}=0$ and $\theta\neq 0$ niederle-morigi-rieger-2016 . We further note that in the SS phase the two SF order parameters $\varphi_{e}$ and $\varphi_{o}$ take different non-vanishing values.

IV.1 Ground-state phase diagram for varying density

Figure 2 shows the ground-state phase diagram as a function of $U_{\infty},\mu,t$, the different colors identify a different phase, the SF phase is the corresponding empty space. In the plane at $U_{\infty}=0$ we recover the mean-field phase diagram of the two-dimensional Bose-Hubbard model fisher-et-al-1989 ; zwerger-2003 . For $0<U_{\infty}/U_{0}<0.5$ the MI lobes shrink along the $\mu$ axis and are sandwiched by CDW phases, which become increasingly visible. Here, the CDW phases are characterized by minimal population imbalance $\theta=1$, corresponding to $n_{e}=n$ and $n_{o}=n+1$ where $n$ is integer, or vice versa. The red regions at the tip of each CDW lobe is SS, the parameter region where the SS phase is different from zero increases with $U_{\infty}$.

Inspecting Fig. 2 we observe that the MI phases vanish at $U_{\infty}=0.5U_{0}$ also for finite tunneling. Moreover, there is a discontinuity at $U_{\infty}=0.5U_{0}$: the CDW phases with population $\{n,n+1\}$ completely disappear and are replaced by CDW phases with maximal population imbalance $\{0,2n+1\}$. This result is in agreement with our analysis in the atomic limit, Sec. II. Moreover, at $U_{\infty}=0.5U_{0}$ and at finite tunneling rate one observes a discontinuous transition from SF to CDW. For $U_{\infty}>0.5U_{0}$ the CDW phases are separated from the SF phase by lattice SS phases, which almost completely surround the tip of CDW regions.

We now consider the values $U_{\infty}=0.3U_{0}$ and $U_{\infty}=0.6U_{0}$ and show the behavior of SF order parameter and $\theta$, respectively, in Fig. 3. We first discuss the case $U_{\infty}=0.3U_{0}$, namely, when the strength of the long-range interaction is below the threshold value $U_{\infty}=U_{0}/2$. In this case, the MI phases are stable. The transition between MI-SF and CDW-SS are characterized by a continuous change of the SF order parameter. However, there is no direct transition between the MI and SS phases. Our numerical results, moreover, predict a direct transition between CDW and MI at $t>0$. The population imbalance changes discontinuously across the CDW-MI transition boundary. In particular, in the vicinity of the transitions between each two insulating lobes, we find a range of parameters where they are metastable: The transition line here corresponds to the parameters where the two states have the same energy. This prediction agrees with the ones of Refs. dogra-et-al-2016 ; sundar-mueller-2016 ; panas-kauch-byczuk-2017 . A direct CDW-MI transition is also predicted by a mean-field treatment in a canonical ensemble wald-et-al-2018 .

We remark that a direct CDW-MI transition, a direct CDW-SF transition, and SS phases at the tip of the CDW lobes have also been found in mean-field studies based on cluster analysis niederle-morigi-rieger-2016 . A further quantitative comparison with the phase diagram reported there is not possible. In fact, the effective strength of the long-range interaction term is not constant across the phase diagram of Ref. niederle-morigi-rieger-2016 , since this depends on the overlap integral between the cavity standing wave and the Wannier functions. There, the Wannier functions are calculated by changing the depth of the confining optical lattice, after which the integrals giving $t$, $U_{0}$, and $U_{\infty}$ are determined.

We now discuss the phase diagrams in Fig. 3 for $U_{\infty}=0.6U_{0}$. Comparison with the left panels show that now the CDW lobes have moved towards smaller chemical potentials, their width (with respect to the chemical potential) has decreased and the critical tunneling rate has increased. The form of this phase diagram qualitatively agrees with the one reported in dogra-et-al-2016 ; sundar-mueller-2016 . Also in this case, discrepancies between the numerics and the analytical lines are visible at the direct transition CDW-SF and at the transition between CDW and SS with different values of the population imbalance. At these phase boundaries the population imbalance varies discontinuously. We discuss in the next section the nature of the other transitions.

Finally, we observe that at fixed chemical potential and at fixed values of $U_{\infty}$ above $U_{0}/2$ the CDW phase has constant population imbalance. This contrasts with the prediction of Ref. liao-et-al-2018 , where a transition between CDW phases with maximal population imbalance as a function of the tunneling rate was reported. Figure 4 shows the occupations $\rho_{e}$ and $\rho_{o}$ of the even and odd sites, respectively, as well as the corresponding SF order parameters as a function of the tunneling rate for the same parameters as in Fig. 6 of Ref. liao-et-al-2018 . We find that in the incompressible phase the population imbalance is constant and equal to $|\theta|=5$. We note that our self-consistent analysis at $t=0$ gives that the CDW with occupations $\{0,4\}$ is metastable with energy $-1.8U_{0}$, while the CDW with occupations $\{0,5\}$ is the ground state with energy $-1.85U_{0}$. Since the mean-field energy does not depend on the tunneling rate in the incompressible phase, then the $\{0,5\}$ CDW is the ground state for all values of $t$ where it is stable. We conclude that a transition like that reported in reference liao-et-al-2018 is not consistent within the static mean field assumption.

Before concluding this section, we briefly compare the phase diagrams in Fig. 3 for $U_{\infty}=0.3U_{0}$ with the ones for dipolar gases, interacting repulsively in two dimensions. Here, mean-field treatments and quantum Monte Carlo calculations report the same phases as for all-to-all coupling, however the ground state phase diagrams are qualitatively different. An important difference is that for dipolar gases there is no direct transition CDW-MI kovrizhin-et-al-2005 ; menotti-trefzger-lewenstein-2007 ; iskin-2011 ; ohgoe-et-al-2012 ; batrouni-scalettar-2000 .

IV.2 Ground state for fixed densities

We now discuss the phase diagram as a function of $t/U_{0}$ and $U_{\infty}/U_{0}$ at fixed density $\rho$. Within our grand-canonical model this implies to find the values of the chemical potential $\mu$, at given $\tilde{t}=t/U_{0}$ and $\tilde{U}_{\infty}=U_{\infty}/U_{0}$, which satisfy the equation $\rho(\mu/U_{0},\tilde{t},\tilde{U}_{\infty})=$ constant. Since the compressibility shall fulfil $\partial\rho/\partial\mu\geq 0$, we can use a bisection algorithm to efficiently find the chemical potential which corresponds to a fixed density. The details are reported in Appendix B. This procedure did not provide a solution for all values of parameters $t/U_{0}$ and $U_{\infty}/U_{0}$ because the compressibility $\partial\rho/\partial\mu$ is not continuous over the full range of $\mu$ values: we find jumps in the density as a function of the chemical potential, as we discuss in what follows.

Figure 5 shows the phase diagram for $\rho=1/2,1,3/2$, and $2$. For $\rho=1/2$ there is no MI phase. Nevertheless, for $U_{\infty}>0$, we observe parameter regions where the ground state is in the CDW phase, corresponding to the occupation $\{0,1\}$ between neighbouring sites. For $U_{\infty}\lesssim 0.1U_{0}$ CDW and SF are separated by a first order phase transition. This phase boundary is characterized by a discontinuity of the population imbalance, the transition line is at a value of the tunneling rate which scales seemingly linearly with $U_{\infty}/U_{0}$ and ends at a tricritical point. After this point the SS phase separates the CDW from the SF phase and the order parameters vary continuously at the transition lines separating SF-SS and SS-CDW. The area enclosed by the dotted lines in the diagram is the parameter region where we could not find any data point, namely, where there is no value of the chemical potential corresponding to $\rho=1/2$. We denote this region by Phase Separation (PS), after observing that simulations for these values in a canonical ensemble using QMC reported negative compressibility and have been linked to a phase separation between the SF and SS phases flottat-et-al-2017 .

We first notice that this phase diagram coincides with the one reported in Ref. dogra-et-al-2016 , apart from the fact that the authors seem to always find a SS phase separating the CDW and the SF phases, and thus they report neither a direct SF-CDW transition nor a PS region. In particular, all transitions they find for $\rho=1/2$ are of second order. This difference, and especially the absence of the PS region, might be attributed to different methods for determining the ground state at a fixed density in a grand-canonical ensemble. The authors of Ref. dogra-et-al-2016 first identify the states at the target density for given $t,U_{\infty}$, and then search for the lowest-energy state in this set dogra-et-al-2016 ; dogra-2017:private-communication . In our work, instead, we first determine the states at the lowest energy as a function of $\mu$ for given $t,U_{\infty}$. In this set of states we then search for the one corresponding to the target density. The PS region corresponds to the parameters for which the target density cannot be found. Details of our analysis are reported in Appendix B.

Remarkably, the plot for $\rho=1/2$ reproduces qualitatively the corresponding diagram obtained with QMC in Ref. flottat-et-al-2017 . In particular, the authors claim to find a direct transition CDW-SF at smaller values of $t/U_{0}$ (larger values of $U_{0}/t$), however they cannot determine its nature due to the fact that the QMC simulations are not conclusive in this parameter regime. The salient difference with our result is that the authors do not report stable SS phases for $\rho=1/2$. This does not exclude, in our view, that a stable SS phase could exist in a small parameter region close to the tricritical point, which might have not been included in the data sampling.

The phase diagrams for $\rho=1$ and $\rho=2$ have a similar structure. For both cases the phases are MI, SF, SS, and CDW with maximal population imbalance. The MI-SF and the SS-CDW transitions are continuous. The MI-CDW, instead, is a discontinuous transition. Moreover, for both densities $\rho=1,2$ there is a direct, discontinuous transition CDW-SF at $U_{\infty}\sim 0.5U_{0}$ which ends at a tricritical point. As $U_{\infty}/U_{0}$ is further increased, this transition line splits into two phase boundaries: the SS-CDW and the SF-SS. The SF-SS transition is continuous except for a small region close to the tricritical point. This region corresponds to the parameter regime for which we find no solution of the equation $\rho(\mu/U_{0},\tilde{t},\tilde{U}_{\infty})=1$. In the case of $\rho=2$, instead, the SF-SS transition is discontinuous close to the tricritical point. On the other hand, the SS-CDW is a continuous transition until a critical value $U_{\infty}(t,\rho)$, after which we find a PS region.

The diagram for $\rho=1$ in Fig. 5 is in full agreement with the one reported in Ref. dogra-et-al-2016 , within the parameter intervals considered. Moreover, it also agrees qualitatively with the phase diagram evaluated using QMC flottat-et-al-2017 , apart for two salient features: The authors do not report a PS and the transition line SS-SF is continuous along the whole branch of their phase diagram. We note that the phase diagram at $\rho=1$ in Fig. 5 is similar to the one of Ref. wald-et-al-2018 , obtained by minimizing the mean-field free energy of a canonical ensemble in a constrained Hilbert space. According to the free energy landscape of Ref. wald-et-al-2018 , the MI-CDW transition is characterized by a large parameter region where the two phases are metastable.

The phase diagram at $\rho=3/2$ in Fig. 5 exhibits a CDW phase with $\{1,2\}$, separated by a discontinuous transition to the CDW phase with $\{3,0\}$ at $U_{\infty}=U_{0}/2$. The CDW $\{1,2\}$ emerges at infinitesimally small tunneling parameters, it has a first-order transition to a SF until a finite value $U_{\infty}<U_{0}/2.$ This CDW$\{1,2\}$-SF phase boundary ends at a tricritical point, after which SF and CDW are separated by a SS phase. The transition SF-SS is continuous in the whole parameter range. On the other hand, SS-CDW$\{3,0\}$ transition, becomes discontinuous for a small interval of values about $U_{\infty}\sim 0.5U_{0}$. Within the SS phase, moreover, there is a discontinuous transition at $U_{\infty}=U_{0}/2$ where the population imbalance undergoes a jump from $|\theta|\approx 1$ to $|\theta|\approx 3$. This jump was reported also in Ref. dogra-et-al-2016 . Instead, QMC studies found it to be a crossover flottat-et-al-2017 . Moreover the direct transition between SF and CDW $\{1,2\}$ seems to not have been found by static mean field calculations in Ref. dogra-et-al-2016 . QMC simulations flottat-et-al-2017 here reported this direct transition, however they could not determine its order. Finally, we notice a region for strong long-range interaction and large tunneling where no solution exists and which was not reported by static mean-field calculations dogra-et-al-2016 .

We note that the phase diagram for $\rho=1/2$ is similar to that of the extended Bose-Hubbard model with repulsive nearest neighbour interaction: For nearest-neighbour interactions and small tunneling rates Quantum Monte Carlo simulations report no stable SS phase but a direct transition CDW-SF batrouni-et-al-1995 ; batrouni-scalettar-2000 ; sengupta-et-al-2005 , which is found to be of first order batrouni-et-al-1995 ; batrouni-scalettar-2000 . Furthermore in the intermediate tunneling regime a SF-SS and SS-CDW transition is observed kimura-2011 ; ohgoe-et-al-2012 . For $\rho=1$, in the extended Bose-Hubbard model with repulsive nearest neighbour interaction and small tunneling rates there is a MI-CDW transition sengupta-et-al-2005 ; kimura-2011 , while the author of ref. kimura-2011 finds a SF-CDW transition for an intermediate tunneling rate $zt=0.3U_{0}$, and SF-SS and SS-CDW transitions for a very large tunneling rate $zt=U_{0}$. However, no PS at $\rho=1$ and $\rho=1/2$ is reported in Refs. kimura-2011 ; ohgoe-et-al-2012 .

V Conclusions

We have performed a mean-field analysis of the phase diagram of the extended Bose-Hubbard model, where the bosons have a repulsive contact interaction and experience an infinitely long-range two-body potential. The systematic comparison between the phase diagram obtained for the cavity Bose-Hubbard model and the one for repulsively interacting dipolar gases in two dimensions shows clear differences already within the mean-field treatment, such as for instance the direct first-order transition CDW-MI at a critical value of the chemical potential, that is absent for the dipolar case. The ground-state phase diagram we calculate mostly agrees with the static mean-field diagram of Ref. dogra-et-al-2016 . There are two important differences: differing from Ref. dogra-et-al-2016 , for the density $\rho=1/2$ and $\rho=3/2$ we predict a direct transition between Superfluid (SF) and Charge-Density Wave (CDW), which is first order for the densities we considered. Moreover, in the region where the authors of Ref. dogra-et-al-2016 predict stable lattice supersolid (SS) phases, we find also regions where instead there is a Phase Separation (PS). We attribute these discrepancies to different methods for determining the ground state at fixed density from the grand-canonical ensemble calculations, as we detailed in Sec. IV.2.

The stability of the SS phase has been extensively analysed by means of Quantum Monte Carlo methods for a canonical and a grand-canonical ensemble in Ref. flottat-et-al-2017 . Our diagrams and the diagrams of Ref. dogra-et-al-2016 at fixed densities, extracted from grand-canonical ensemble calculations, are in remarkable qualitative agreement with the QMC diagrams in the interval of parameters where the QMC diagram have been determined. The discrepancies regard the stability of the SS regions at fixed densities. These discrepancies could be due to the fact that the QMC collected data did not sample the regions where these differences are found (as one could conjecture by taking the parameters reported in Ref. flottat-et-al-2017 for which the stability of the SS phase was analysed and mapping them into our phase diagram). Most probably, the discrepancy arises from the fact that our static mean-field approach cannot appropriately take into account the interplay between strong long-range interactions and quantum fluctuations.

Acknowledgements.

The authors are thankful to George Batrouni, Nishant Dogra, Rosario Fazio, Gabriel Landi, Astrid Niederle, Heiko Rieger, Katharina Rojan, and Sascha Wald, for discussions. Financial support by the DFG priority program no. 1929 GiRyd, by the DFG DACH program ”Quantum crystals of matter and light”, and by the German Ministry of Education and Research (BMBF) via the Quantera project “NAQUAS”. Project NAQUAS has received funding from the QuantERA ERA-NET Cofund in Quantum Technologies implemented within the European Union’s Horizon 2020 Programme. CC acknowledges funding from the Alexander-von-Humboldt Foundation for her visit to Saarland University and grant BID-PICT 2015-2236.

Appendix A Numerical calculation of the ground state

In this Appendix, we describe the algorithm used to find the self-consistent ground state of the local mean field Hamiltonians (17). We measure all physical parameters of the Hamiltonian in units of the on-site interaction, $\tilde{\mu}=\mu/U_{0}$, $\tilde{U}_{\infty}=U_{\infty}/U_{0}$, $\tilde{t}=zt/U_{0}$ and obtain the Hamiltonians

[TABLE]

with the same eigenenergies and -states as the Hamiltonians (17). We fix the parameters of the Hamiltonian $\tilde{t}$, $\tilde{U}_{\infty}$ and $\tilde{\mu}$. The mean-field order parameters $\varphi_{\text{e}}$, $\varphi_{\text{o}}$ and $\theta$ are now the free variables. The problem is formulated as follows. We first introduce the function

[TABLE]

where $\langle\cdot\rangle_{s}$ denotes the single-site expectation value with respect to the ground state of the Hamiltonian $H_{s}(\varphi_{\text{e}},\varphi_{\text{o}},\theta)$, for $s\in\{\text{e},\text{o}\}$. Further, we define $F$ to be the set of fixed points of $f$,

[TABLE]

The goal is to find the self-consistent order parameters which minimize the energy per site,

[TABLE]

The basic idea of the algorithm is that of fixed point iteration: Apply $f$ repeatedly to some random $(\varphi_{\text{e}},\varphi_{\text{o}},\theta)$, until applying it again does not significantly change the input habibian-et-al-2013 ; dhar-et-al-2011 ; niederle-morigi-rieger-2016 .

We measure the distance between mean-field order parameters by the infinity-norm and relax the criterion for $(\varphi_{\text{e}},\varphi_{\text{o}},\theta)$ to be a fixed-point to

[TABLE]

where $(\varphi_{\text{e}}^{\prime},\varphi_{\text{o}}^{\prime},\theta^{\prime})=f(\varphi_{\text{e}},\varphi_{\text{o}},\theta)$, and $\varepsilon$ is some predefined tolerance, e. g. $\varepsilon=10^{-6}$.

This naïve algorithm has the following problems, however: First, if the algorithm converges to some point, there is no guarantee that this point minimizes the energy per site. Second, the algorithm is not guaranteed to converge. Third, the algorithm sometimes converges sublinearly, and thus extremely slowly.

We approach the first problem by taking a sufficient large number of initial guesses.We always deterministically take the following $75$ initial guesses: $(\varphi_{\text{e}},\varphi_{\text{o}},\theta)\in\cup_{\{n\in 0,\ldots,24\}}\{(0,0,n),(0.001,0.002,n),(0.1,0.2,n)\}$. Additionally, we use the Mersenne Twister and Ranlux48 algorithms to pseudorandomly sample initial guesses from the $\mathrm{Cauchy}(0,1)$ distribution. We find that $50$ random initial values are sufficient to find the minimal energy and verify this by taking more initial guesses and verifying that the energy does not decrease significantly.

The problem that the algorithm sometimes does not converge manifests in the way of cycles of the form that $f(\varphi_{\text{e}},\varphi_{\text{o}},\theta)\approx(-\varphi_{\text{e}},-\varphi_{\text{o}},-\theta)$ and $f(-\varphi_{\text{e}},-\varphi_{\text{o}},-\theta)\approx(\varphi_{\text{e}},\varphi_{\text{o}},\theta)$. We detect this by comparing not only the mean-field order parameters, but also their absolute values. If the difference of the absolute values is smaller than $\varepsilon/10$ for $1000$ consecutive iterations, we re-run the algorithm with the absolute values of the final order parameters $(\varphi_{\text{e}}^{\prime},\varphi_{\text{o}}^{\prime},\theta^{\prime})$ as an initial guess. To ensure that we still find the minimal energy, we compare the energies of the result of the initial run $(\varphi_{\text{e}}^{\prime},\varphi_{\text{o}}^{\prime},\theta^{\prime})$, with that of the second run, $(\varphi_{\text{e}}^{\prime\prime},\varphi_{\text{o}}^{\prime\prime},\theta^{\prime\prime})$. For this comparison, we do not consider the eigenvalues of the two Hamiltonians, but the expectation value of the updated Hamiltonian with respect to the ground state of the Hamiltonian before the update. If the energy of the second run is smaller, we accept this solution, otherwise we reject it.

Finally, we note that the algorithm converges to the set tolerance within a few hundred or thousand iterations (i. e. applications of $f$) in a large region of the phase diagram. In this case, the algorithm converges linearly. However, in some cases it converges sublinearly and extremely slowly. Figure 6 shows a comparison of two cases, for two points in the phase diagram which are close to each other, and identical initial guesses.

The algorithm converges that slow only at relatively few points in the phase diagram. We have verified that the number of iterations does not strongly influence the value of the resulting mean-field order parameters, by comparing the results after $10^{4}$ and $10^{6}$ iterations.

For finding the ground state of the Hamiltonians (46), we truncate the Hilbert space of each site taking the cutoff $n_{\text{max}}=31$, leaving us with two tridiagonal real symmetric $32\times 32$ matrices (in the Fock basis), which we diagonalize numerically. We identify the cutoff $n_{\text{max}}=31$ by performing calculations also for $n_{\text{max}}=23$ and $n_{\text{max}}=63$ and verifying that the results do not differ significantly.

We implemented the algorithm in the C++ language, using the GNU compiler collection (versions 7.2 and 8) and clang (versions 5 and 6). For diagonalization, we used the library Eigen3. We verified the self-consistency of the results with a partial implementation of the algorithm in the Python language (version 3.6) and using the numpy library (version 1.14.5) for diagonalization. software

Appendix B Supersolidity and phase separation for fixed densities

In this Appendix, we report details of the calculations for determining the phase diagram for constant densities of Fig. 5.

We obtain the order parameters for fixed densities by adjusting the chemical potential such that the ground state has the target density. More precisely, we perform a bisection algorithm starting at $\mu=-U_{0}$ which gives a density lower than the target density and $\mu=3U_{0}$ which gives a density higher than the target density. In every step, the ground state for the midpoint of the $\mu$ interval is calculated following the procedure detailed in Appendix A, its density is computed, and the interval is halved such that at the lower (upper) point of the interval, the density is smaller (larger) than the target density. We repeat this up to twenty times until either the target density is reached or we conclude that a solution is not possible. When the density is not attained with the required precision, which we set to $\varepsilon_{\rho}=10^{-4}$, we name the corresponding point of the phase diagram phase separation, following Ref. flottat-et-al-2017 . Otherwise, we determine the phase from the values of the order parameters $\theta$ and $\varphi$.

For a density of $\rho=0.5$, we find both a SS region and a region of phase separation, unlike Refs. dogra-et-al-2016 ; flottat-et-al-2017 . Fig. 7 shows the $\rho(\mu)$ curve for a PS point. Fig. 8 shows the $\rho(\mu)$ curve for a different point, where the density $\rho(\mu\approx-0.132)=0.5$. The same figure shows the superfluid order parameter and the even-odd imbalance; for the parameters where $\rho(\mu)=0.5$, the ground state is SS.

In Fig. 9, we show the $\rho(\mu)$ curve for a PS point for density $\rho=1$, which was not reported so far dogra-et-al-2016 ; flottat-et-al-2017 . As shown in fig. 5, we also find a PS region for $\rho=1.5$, which was not reported so far dogra-et-al-2016 ; flottat-et-al-2017 .

In Fig. 10, we show the order parameters along cuts of the phase diagram in figure 5, specifically the density $\rho=0.5$. We show similar plots for the density $\rho=1$ in figure 11. We used plots similar to the ones shown in Figs. 10 and 11 to determine all phase boundaries of Fig. 5.