Light-induced quantum droplet phases of lattice bosons in multimode cavities

Multimode optical cavities can be used to implement interatomic interactions which are highly tunable in strength and range. For bosonic atoms trapped in an optical lattice, cavity-mediated interactions compete with the short-range interatomic repulsion, which we study using an extended Bose-Hubbard model. Already in a single-mode cavity, where the corresponding interaction has an infinite range, a rich phase diagram has been experimentally observed, featuring density-wave and supersolid self-organized phases in addition to the usual superfluid and Mott insulator. Here we show that, for any finite range of the cavity-mediated interaction, quantum self-bound droplets dominate the ground state phase diagram. Their size and in turn density is not externally fixed but rather emerges from the competition between local repulsion and finite-range attraction. Therefore, the phase diagram becomes very rich, featuring both compressible superfluid/supersolid as well as incompressible Mott and density-wave droplets. Additionally, we observe droplets with a compressible core and incompressible outer shells.

Introduction.Ultracold atoms in optical cavities represent a unique experimental platform for the study of strongly interacting quantum many-body systems of light and matter [1,2].The cavity-mediated interactions are naturally long-ranged and can easily overcome the shortrange interatomic repulsion.The competition between these two types of interaction can be enhanced by trapping the atoms in an optical lattice.Already for a single-mode cavity, which mediates global-range interactions, a rich phase diagram has been experimentally observed, featuring, besides the known superfluid and Mott-insulating phases, a lattice supersolid as well as an incompressible density-wave phase [3].Such experiments implement a Bose-Hubbard (BH) model extended by a global-range sign-changing interaction, which has been intensively studied in recent years [4][5][6][7][8][9][10][11][12][13][14].
More recently, the coupling to multiple modes of the cavity has emerged as a very promising route towards an even richer many-body physics.Already the use of two cavity modes has allowed to observe a supersolid with continuous translational symmetry breaking [15][16][17][18].By tuning cavitites around the confocal degenerate point [19,20], it has recently become possible to even control the range of the cavity-mediated interaction [21,22], and to realize an optical lattice featuring phonons [23].With finite-range cavity-mediated interactions, supersolids have been predicted to feature crystalline topological defects and to appear through nonmean-field phase transition dominated by fluctuations [24,25].
Here we demonstrate that, due to the competition between cavity-mediated finite-range attraction and the onsite repulsion, the phase diagram is actually dominated by a variety of droplet phases, i.e. self-bound manyparticle quantum objects, while the extended phases exist only for sufficient repulsion (see [26] for the classical, purely attractive case).We study a BH model ex-tended by an attractive tunable-range interaction, with and without additional sing-changing modulation, which can be realized in state-of-the-art confocal-cavity experiments.We compute the ground-state phase diagram with a worm quantum Monte Carlo algorithm, using a canonical version [27][28][29] necessary to describe droplet phases.We observe a complex sequence of transitions between droplets of different sizes, and of compressible (superfluid or supersolid) as well as incompressible (Mott or density-wave insulating) nature, governed by the competition between the local repulsion and the finite-range attraction.Within the superfluid-droplet phase, we additionally observe a sequence of crossovers between the fully superfluid droplets and droplets with superfluid core and Mott-like outer shells.
Lattice bosons with dipolar interactions [41,42] and mixtures [43,44] have been proposed for implementing various types of extended Hubbard models [45].While currently for dipolar gases and bosonic mixtures the challenge is to reach appreciable interactions beyond the nearest neighbour, cavity-mediated interactions are instead naturally strong at large distances.This feature is shared by other types of photon-mediated interactions based on refraction [46] or diffraction [47][48][49] (the latter also predicted to induce droplet phases [50,51]), which however lack the tunability of the range.
Model.We consider ultracold bosons trapped in an optical lattice in the two different configurations illustrated in Fig. 1.A quasi-1D geometry with monotonous cavitymediated interactions, and a quasi-2D with sign-changing FIG. 1.We consider ultracold bosons in two configurations corresponding to two different cavity-mediated interaction potentials.The latter are generated in the dispersive regime via two-photon transitions involving retroreflected lasers beams (red arrows) which are red detuned from a given degenerate family of a multimode cavity.The interaction potential results from the interference between the laser and the cavity field.In (a), a quasi 1D gas is trapped in an optical lattice within the midplane of the cavity: z = 0, along the x direction perpendicular to the laser.This results into the monotonous interaction depicted in (c), which decays over a scale ξ controlled by the number of degenerate, transverse cavity modes [21,52].In (b), the lattice gas is quasi-2D and two perpendicular retroreflected laser modes interfere with the cavity field.This induces the sign-changing interaction depicted in (d).The trapping optical-lattice potentials are generated by three additional far-off detuned lasers (oriented along the x, y, and z directions, not shown) which do not interfere with the cavity field [3].
interactions.As derived in detail in the Supplementary Material, the system can be described by an extended BH model Here i, j denotes nearest-neighbor sites in a D dimensional square lattice, t is the hopping between the nearest neighboring sites of the optical lattice, b † i and b i are bosonic creation and annihilation operators, n i = b † i b i , U > 0 is the on-site contact repulsion, and V i, j is the tunable-range interaction mediated by the cavity between particles at sites i and j.We consider a multimode optical cavity, where the finite range ξ is achieved by having a large number of transverse modes within a given degenerate family.Alternatively, Floquet modulation [53] or multifrequency driving [54] can render the range finite even for non-degenerate cavities.Experiments with ultracold bosons in confocal cavities have demonstrated the ability to tune this range [21].We model this using the following generic form of the interaction potential: where V 0 > 0 (locally attractive interaction).The signfactor s is either s i = 1 for the monotonous attractive interaction or s i = (−1) (ix+iy) for the sign-changing interaction.All lengths are measured in the unit of lattice spacing.We stress that the results presented in the following do not qualitatively depend on the particular form of the envelope decay, as long as a single scale ξ can be associated with it.
The length scale ξ introduced above determines the characteristic radius of the droplet and helps to define a proper thermodynamic limit for a droplet phase, which is given by L → ∞, N → ∞, ξ → ∞, V 0 → 0, while V 0 ξ D = const and N/ξ D = const n 0 .Here V 0 ξ D gives the characteristic integrated interaction strength and the "droplet filling" n 0 gives the characteristic density of particles in a droplet of radius ξ.A further important parameter is the characteristic slope of the density profile, σ dn/dr n 0 /ξ = N/ξ D+1 (see Supplementary Material).This steepness parameter vanishes in the droplet thermodynamic limit defined above, and quantifies the amount of finite-size effects in the phase diagram.As we discuss later and in the Supplementary Material, these finite size effects are experimentally measurable and it is important to understand them for the correct extrapolation to the thermodynamic limit.
Method.We map the ground-state phase diagram of the model using a canonical worm Quantum Monte Carlo algorithm [27][28][29] with fixed total number of particles N = i n i and periodic boundary conditions.Since we explicitly work in the canonical ensemble, we can deal with arbitrary occupation numbers (up to the total number of particles in the system) and thus for the search of potential droplet phases the method is superior with respect to other traditional methods of studying the BH model, like grand-canonical quantum Monte Carlo [55][56][57] methods or DMRG [58], which require fine tuning of the chemical potential and/or cut-offs for the local occupation.In order to construct the ground-state phase diagram of the model, we perform calculations at small non-zero temperature T = 0.1t (and check that this temperature is low enough by comparing our method with Lanczos exact diagonalization up to L = N = 16 systems in 1D).
We characterize the droplets using two observables: the droplet order parameter O dr = where N s = L D , and the single-particle superfluid correlations between neighboring sites b † i+1 b i (in 1D monotonous interaction case) and next-neighboring sites b † ix+1,iy+1 b ix,iy (in 2D sign-changing interaction case), averaged over all sites i. Monotonous interaction.Figure 2 shows the groundstate phase diagram for the case of monotonous interaction s i,j = 1 in 1D (setting of Fig. 1(a)).

SF-dr
For U = V 0 = 0 the system is in the SF state.With increasing V 0 we enter a single-site-droplet phase at V 0 /t ∼ 1/N [59].Upon increasing U and decreasing V 0 we observe smooth crossovers between droplets of different "diameters" (Fig. 2(a)), up to the maximum diameter bounded either by the size of the system or the total number of particles d max min(L/2, N ).Droplets whose diameter is larger than ξ are only weakly bound as the outer particles are not bound to the center.After the maximum droplet size is reached, a first-order transition to a phase with uniform density occurs, either a Mott insulator for integer fillings or a SF for non-integer fillings.
Fig. 2(b), shows the nearest-neighbor superfluid correlators, where we observe the various superfluid and Mott quantum droplet phases.The Mott phases here correspond to different fillings at the core of the droplet, e.g.
The existence of the transitions between SF and Mott droplets can be qualitatively understood by means of the standard Bose-Hubbard phase diagram [60], whereby the chemical potential is set by the local occupation at the droplet core.For V 0 = 0 and a fixed integer average filling ν = N/L, upon increasing U we have a conventional SF-Mott transition at U/t ≈ 4ν [61].Analogously, for the droplet phases, with increasing U along lines V 0 /U = const the superfluid correlator between neighboring sites shows SF-Mott transitions.We emphasize however that, unlike the case of the Mott-SF transitions in an external confining potential, here the droplet size and density distribution is not externally fixed but rather emerges from the competition between local repulsion and finite-range attraction.Moreover, the droplet spontaneously breaks the translation symmetry of the Hamiltonian.
Additionally, we compute the real-space density profiles of the droplets (insets 1 ○-7 ○ of Fig. 2(b)).Moving outwards in the radial direction of the phase diagram we first observe the fully superfluid droplet 1 ○ (see Supplementary Materials for an analytic Thomas-Fermi description), which gradually develops Mott-type outer shells, and finally crosses over to a droplet with only the core being superfluid 4 ○. Moving along the vertical line U/t = 50 we observe instead transitions between SFdroplets 4 ○, 6 ○ and Mott droplets 5 ○, 7 ○.
The inset 8 ○ of Fig. 2(b) shows a zoom in the V 0 /t < 0.5 part of the phase diagram, featuring the extended SF phase.Using (1) we can estimate the critical value of V 0 for the transition between the Mott droplet with n c = 1 and the SF-phase.In the former phase the energy is dominated by the interaction part E Mott−dr ∼ −N V 0 ξ D , while in the latter phase the energy is dominated by the kinetic part E SF ∼ −N t.Equating the two energies, the transition value can be estimated as V 0 ∼ t/ξ D .
We have explicitly checked (finite-size scaling data from L = 100, N = 20, ξ = 25 to L = 500, N = 100, ξ = 125 is not shown) that the droplet regime of the phase diagram has a well-defined thermodynamic limit, L, N, ξ → ∞, while N/ξ = const.For finite-size droplets featuring a steepness parameter σ 1, we observe additional features, e.g.mesoscopic first-order phase transitions associated with change in the droplet size (see Supplementary Material).Similar examples are visible in the 2D calculations shown next.
Sign-changing interaction.Now we turn to the case of the sign-changing interaction, which makes possible the existence of the density-wave (DW) and supersolid (SS) phases.We consider the two-dimensional setting of Fig. 1(b).
The infinite-range case, ξ = ∞, has been studied experimentally in a single-mode cavity [3].The phase diagram we find in the case of a finite range, shown in Fig. 3 for the relatively short range ξ = 1 and unit filling n = N/L 2 = 1, features the same extended phases as the infinite-range case (studied with QMC in [9]), plus quantum droplet phases, which occupy a substantial part of the phase diagram.Note that in order to perform the comparison between the ξ = ∞ and finite ξ cases we use the rescaled coordinate of vertical axis, V 0 ξ D /t.
For V 0 = 0 the system is either in the SF phase (at U/t 15) or the Mott phase (at U/t 15, Fig. 3).Increasing V 0 favors the density-modulated SS and DW phases.In these phases we can estimate an effective hopping within a given sublattice using second-order perturbation theory.In extended SS and DW phases, as well as for droplets with d ξ, the effective hopping between sites of the even (odd) checkerboard sublattice is t eff ∼ n 0 t 2 /(γn 0 V 0 ξ D −(n 0 −1)U ), with γ some numerical constant.For smaller droplets with d ξ the corresponding estimate gives t eff ∼ n 0 t 2 /(γN V 0 −(n 0 −1)U ).In both cases, upon moving in the radial direction in the phase diagram the extended SS and the SS-droplet phases feature much smaller SF correlators in comparison to the case of monotonous interactions, where b † i+1 b i ∼ t/U .This is due to the fact that for the density-modulated phases b † ix+1,iy+1 b ix,iy ∼ t eff /U ∼ t 2 /U 2 acquires an additional t/U suppression factor for increasing U, V 0 ξ D and U/V 0 ξ D = const.Hence, the lobe structure for different DW-droplet phases is not visible in Fig. 3(a).Still, t eff is finite in the thermodynamic limit N, ξ → ∞, V 0 → 0, since n 0 = N/ξ D , V 0 ξ D , as well as their product, N V 0 , are constant.
In order to make the lobe structure of the quantum droplet phases visible, in Fig. 4 we show the phase diagram for a small lattice filling n ≈ 0.08 and low steepness parameter σ ≈ 0.16 (see Supplementary Material for the finite-size case with large σ).Like for the monotonous interaction, we observe a non-trivial change in the superfluid correlator in the azimuthal direction (Fig. 4(b)), due to appearance of the self-consistent DW-droplet lobes separated by regions of superfluid droplets.The insets of Fig. 4(c) show typical Mott-droplet real-space density profiles.We note that here σ is not small enough to eliminate finite-size effects: in the upper-left corner of Fig. 4(a) we see a sequence of first-order transitions between superfluid droplets of different sizes.
Conclusions and outlook.In this work, we have proposed an experimental setting for the creation self-bound quantum droplets for bosonic particles in an optical lattice, where the finite-range attractive tunable interaction is induced by a multimode optical cavity and competes with the local repulsion.We modelled such a system by a rather generic class of extended Bose-Hubbard models.We observed various types of lattice droplet phases: superfluid/Mott droplets for the monotonous finite-range interactions and supersolid/density-wave droplets for sign-changing interactions.Unlike the case of lattice bosons in a trap, droplets are self-bound objects which break the translation symmetry spontaneously, so that the local chemical potential emerges from the competition between local repulsion and finite-range attraction.
Other unexplored classes of extended BH-models related to the one studied here might be also realized by trapping dipolar Bose gases [30,31,[34][35][36][37] or bosonic mixtures [32,33] in a lattice.The numerical methods and the analysis employed here should offer a helpful guide for the future investigation of new quantum droplet phases in such models.
Under realistic conditions for large local densities of particles also three-body effects, not considered here, can play an important role and show up either as a three-body interactions [45,62,63], favoring extended droplets, or three-body losses, making the system metastable, like in the recent dipolar gas experiments [35][36][37].
Finally, we mention that the droplet physics discussed here is rather generic and our analysis might also be relevant in other contexts, for example for the quantum simulation of black holes as real-space Bose-Einstein condensates of gravitons [64].In this Supplementary Section we perform the Thomas-Fermi analysis of the superfluid droplets and find their density profiles.
Here we will work in the continuum limit: i.e. the limit when all characteristic length scales of the problem such as the range of the interaction potential and the radius of the droplet ξ, R 1 (let us recall that all lengths are measured in the units of lattice constant a), and the local occupation numbers n 1 in the core of the the droplet.For simplicity, we will work in one dimension.
The density profile n TF ≡ n(x) in the Thomas-Fermi approximation in a self-consistent form is given by where )dx is the effective potential experienced by a particle placed at x, V (x) = −V 0 exp(−x 2 /ξ 2 ) is the interaction potential, and θ is the Heaviside step function.Instead of solving the above integral equation we choose an ansatz for the density profile and optimize it with respect to the parameter R (the radius of the droplet).Normalization is chosen in such a way that R −R n(x)dx = 1.Note that we can systematically improve the ansatz by introducing the higher even powers of x and the corresponding variational parameters.
The total Thomas-Fermi energy of the system can be written in the continuous limit as where On the one hand, it gives an estimate for the number of SF-and Mott-droplet stripes that can fit within a droplet phase with fixed size.On the other hand, it gives an estimate of the spatial derivative of the density of the droplet σ ∼ ∂n(x)/∂x ∼ (N/ξ)/ξ = N/ξ 2 .We thus name σ the droplet steepness parameter.According to it we classify the droplets as steep (σ 1) or flat (σ 1).Note that in the thermodynamic limit for the droplet, N, ξ → ∞, with N/ξ = n 0 = const, the steepness parameter σ = n 0 /ξ → 0, so only the flat droplets survive.In this regime eq.(S12) (except of the last proportionality specific to the Thomas-Fermi limit) still gives a good estimate for the characteristic width of SF/Mott droplet regions, from which we see that all the phases with the droplet radius scaling as R ∼ ξ survive in the thermodynamic limit.At the same time, all the regions of the phase diagram corresponding to the droplets of finite fixed size R shrink in the thermodynamic limit as ∆V Mott−SF 0 ∼ 1/N (see first equality in eq. ( S12)).
In the main text we study the only the flat droplet regime relevant to the thermodynamic limit, while the steep droplets might be present in finite systems also relevant for experiments and we discuss them in the next Supplementary Section.

II. FINITE-SIZE REGIME OF STEEP DROPLETS A. Monotonous interaction in 2D
Performing the same analysis done for the flat-droplet regime in the main text, we consider here the steep droplet regime, both for the monotonous interaction in 1D and for the sign-changing interaction in 2D.First we consider the steep droplet regime for the monotonous interaction in 1D.For U = 0, V 0 > 0 we are the single-site droplet phase.Upon increasing U and decreasing V 0 we observe a sequence of first-order phase transitions between droplets of different "diameters" (Fig. S2(a)), up to the maximum diameter bounded either by the size of the system or the total number of particles d max min(L/2, N ) = 15.We observe that, as soon as d exceeds min(ξ, N ), the corresponding region of the phase diagram shrinks with increasing d, indicating the fact that d = min(ξ, N ) is the characteristic diameter of the droplet.
For steep droplets, each region of the phase diagram corresponding to a droplet of a given size further subdivides into smaller regions (Fig. 2(b)).These finer transitions correspond to a redistribution of the droplet density profile while its diameter stays fixed.Between two stable configurations with almost non-fluctuating local particle densities (Mott states), which correspond to different discretizations of the optimal continuous droplet profile (upper and lower density profiles in Figs.2(c)), there is always a superfluid region with local density fluctuations.
The steepness parameter σ determines how many different local occupations of the core are compatible with a given droplet size (Fig. 2(b)), that is, how many Mott lobes fit into a given sector of fixed d in Fig. 2(a) (see the discussion in the previous Supplementary Section).Now we turn to the steep droplet regime for the sign-changing interaction in 2D.Like for the monotonous interaction, for each droplet phase we observe a structure of superfluid radial lines in the phase diagram, i.e. a non-trivial change in the superfluid correlator in the azimuthal direction (Fig. S3(b)).It is due to the fact that different, one or another discretization of the density profile n(r) becomes more favorable, giving rise to different stable Mott structures (like in Fig. S3(c)).In the middle of the crossover between two such structures we observe an increase in the superfluid correlators.

III. DERIVATION OF THE EXTENDED BOSE-HUBBARD MODEL A. Cavity-mediated interaction
In this Supplementary Section we describe the setting we propose for the experimental realization of the interaction we use in the main text.Apart from slight variations in the laser-beam orientation, the setup we propose is the one implemented with ultracold bosons in a confocal cavity in Stanford [21][22][23]52].
We consider an ensemble of atoms treated as two-level systems with level splitting ω a , placed in a multimode optical cavity.The atomic transition is off-resonantly driven by either one or two laser beams oriented in the direction transverse to the cavity axis.Each beam has a frequency ω L and the total laser field is indicated by Ω(r).The cavity possesses (nearly) degenerate modes with frequencies ω cα (where α is the mode index).The coupling between the atomic transition and a given cavity mode is given by g α (r) = g 0 Ξ α (r), where Ξ α (r) is the normalized mode function of the cavity mode.We work in the dispersive regime, where the detuning ∆ cα = ω cα − ω L between the laser and the cavity degenerate family is much smaller than the detuning ∆ a = ω a −ω L between the laser and the atomic transition.In this case, we can adiabatically eliminate the atomic excited state and work only with itinerant atoms all the time in the electronic ground state.The optical potential that the atom experience is induced by two-photon transitions [2].The latter can either involve two laser-photons and give rise to a static spatially-dependent Stark shift or involve one photon from the laser and one from the cavity, which induces a dynamical optical potential whose spatial dependence results from the interference between the laser and the cavity field.We neglect the dynamical potential resulting purely from cavity photons as it can be suppressed by increasing the detuning ∆ a [21].Assuming now that the cavity detuning ∆ cα is larger than the dispersive coupling and the atomic motional scales we can further eliminate the cavity field adiabatically and obtain the following cavity-mediated interaction [2,21,22,52] In order to have access to superradiant phase transitions and the corresponding spatial ordering, the laser must be red-detuned from the cavity degenerate family ∆ cα > 0. In this case the interaction is attractive at small distances r → r i.e.V (r, r) < 0, which creates the possibility for droplet-formation.
In the main text, we consider the generic ansatz Eq.( 2) for the interaction potential, decaying with the characteristic length ξ.Such translationally-invariant potentials can be directly created by concentric optical cavities, which are however experimentally not available in combination with an ultracold atomic cloud.Confocal cavities are instead already available in Standford and have been shown to mediated tunable range interactions [21].Differently from our ansatz, the interaction in that case contains an additional non-translationally-invariant contribution of the type cos(rr ).However, using a second laser near resonant with a second nearby degenerate family of the same parity, this term can be eliminated, as has been experimentally demonstrated in [23].Finally, there is a further contribution not included in our model, which comes from the interaction of the atoms with their images on the other side of the cavity axis.This contribution can be also eliminated by placing the cloud sufficiently away (at distances ξ) from the cavity axis.This has been also experimentally demonstrated in [23].

B. Lattice geometry
We consider the situation where the atomic cloud is trapped within a deep optical lattice formed by an additional set of laser beams, building the field Ω ext (r), with a detuning ∆ ext = ω a −ω ext from the atomic transition which is sufficiently different from ∆ a to avoid interference with the lasers involved in the interaction.Such a condition is achieved in state of the art experiments with ultracold bosons in a cavity [3].In what follows, we will assume that |V ext (r)| |V L (r)| and hence neglect the latter.

One-dimensional setting
We choose the z-axis to be parallel to the cavity axis.We will use one laser beam for inducing the cavity-mediated interactions and three lasers for creating the external optical lattice that confines the gas to well-separated 1D chains, parallel to the x-axis (see Fig. 1 in the main text).The external optical potential reads If Ω yz ext Ω x ext , the two latter terms confine the gas to the tubes y = z = 0. a) Sign-changing interaction.The interaction-inducing laser is oriented along the x-axis: Ω(r) = Ω cos(kx).The resulting interaction potential reads Let k ext ≈ k = 2π/λ.For ∆ ext > 0 the external potential confines the gas near points kx = πn ⇔ x = nλ/2 (n ∈ Z).
Near such points r n = (nλ/2, 0, 0) we have , which is the desired 1D sign-changing interaction.b) Monotonous interaction.Now we orient the interaction-inducing laser along the y-axis.Since we consider a deep trapping lattice such the atoms are strongly confined along y, y = 0, the sign-changing part of the interaction is approximately constant: cos(ky) cos(ky ) = 1 and we obtain the 1D monotonous attractive interaction considered in the main text: We chose the polarizations of the interaction-inducing lasers to be e 1 = e y ⊥ e 2 = e x , so they can scatter into the cavity.Let us also first assume that for the optical lattice lasers we have ẽ1 ⊥ ẽ2 ⊥ ẽ3 and Ω z  ), so the interaction potential is translation invariant on the set of all points of the optical lattice as well as all its subsets.We show the periodicity of this interaction in Fig. S4(a): if both particles sit on the same checkerboard sublattice their interaction energy is either positive or negative; if they sit on the different sublattices, their interaction energy is zero.We denote this case as Z 4 interaction.b) Monotonous interaction.If we choose instead the frequencies of the external optical lattice lasers and the interaction-inducing lasers in such a way that ω ext = 2ω (so, k ext = k/2), then the period of the external optical lattice would be doubled and we obtain a translation-invariant monotonous interaction (Fig. S4(b)).c) Sign-changing interaction.This is obtained by choosing ẽ1 ||ẽ 2 ⊥ ẽ3 and k ext ≈ k.In this case thus the lattice lasers 1 and 2 interfere, so that V ext (x, y, 0) = −

8 FIG. 2 .
FIG. 2. The ground-state phase diagram of the 1D system with monotonous interactions in the coordinates V0 (long-range attraction) vs U (on-site repulsion).The droplet thermodynamic limit ξ, N → ∞ with n0 = N/ξ = 0.8 is calculated by extrapolation up to L = 500, N = 100, ξ = 125.(a) Droplet order parameter O dr .(b) Nearest-neighbor superfluid correlator b † i+1 bi .Insets 1 ○-7 ○ show the spatial density profiles of the droplets at the corresponding points in the phase diagram; the error bars show the standard deviation of density from its mean value.Inset 8○ shows the SF-correlator for in the region V0/t < 0.5 featuring the SF phase.

FIG
FIG.S1.The Thomas-Fermi picture for the monotonous interaction in 1D.(a) The Thomas-Fermi energy ETF (in arbitrary units) versus the radius of the droplet R/ξ, given by eqs.(S3)-(S6) for V0ξ/U = 0.6 (droplet state with radius R0) and 0.5 (uniform state).For the actual system with discrete fillings the latter case corresponds to the droplet with local density = 1 particle per site (average density < 1 is assumed).(b) The Thomas-Fermi density profiles n(x) for the for different values of V0ξ/U .Here we use ξ = 25, the density is calculated by the imaginary time propagation of GPE with small hopping t = 0.01V0, and normalization n(x)dx = 1.
FIG. S2.Steep droplet regime for monotonous interaction in 1D.Low-temperature phase diagrams in the coordinates V0 vs U for σ = N/ξ 2 ≈ 33 1 (L = 30, N = 300, ξ = 3).(a) Droplet order parameter O dr .(b) Nearest-neighbor superfluid correlator b † i+1 bi .(c) Spatial profiles of the droplets corresponding to the points in the phase diagrams Fig. (b).The errorbars show the standard deviation of density from its mean value -they are stretched by factor of 100 for the visual clarity.(d) Schematic standard Bose-Hubbard phase diagram.Dashed lines in all figures are guides for the eye for the droplet-size change mesoscopic phase transitions, vanishing in the thermodynamic limit.Solid lines show the transitions present also in the thermodynamic limit.SF phase in (a), (b) are in the region U/t, V0/t 10 (not seen at the given scale).
FIG. S3.Steep droplet regime for sign-changing interaction in 2D.Low-temperature phase diagram in the coordinates V0 vs U for σ = N/ξ 3 = 6.4 (L = 10, N = 100, ξ = 2.5).(a) Droplet order parameter O dr (the insets show spatial structures of the droplets).(b) Superfluid correlator b † ix+1,iy +1 bi x ,iy .(c) Spatial structures of the droplets at the corresponding points of phase diagram; the numbers show the local occupation numbers in the Mott regime (for the unmarked sites n i = 0).

2 a(e 1
FIG.S4.Periodicities of the interaction potential (red and blue circles) and of external optical lattice (black circles) in the 2D setting for (a) the Z4 interaction, (b) the monotonous interaction, and (c) the sign-changing interaction considered in the main text. ext

Ω xy ext 2 ∆a
(cos(k ext x) + cos(k ext y))2 and we obtain the lattice configuration shown in Fig.S4(c).It corresponds to the 2D sign-changing interaction considered in the main text.Note that the optical lattice axes xỹ are 45 • rotated with respect to the xy axes.