Bose-Hubbard model with occupation dependent parameters

We study the ground-state properties of ultracold bosons in an optical lattice in the regime of strong interactions. The system is described by a non-standard Bose-Hubbard model with both occupation-dependent tunneling and on-site interaction. We find that for sufficiently strong coupling the system features a phase-transition from a Mott insulator with one particle per site to a superfluid of spatially extended particle pairs living on top of the Mott background -- instead of the usual transition to a superfluid of single particles/holes. Increasing the interaction further, a superfluid of particle pairs localized on a single site (rather than being extended) on top of the Mott background appears. This happens at the same interaction strength where the Mott-insulator phase with 2 particles per site is destroyed completely by particle-hole fluctuations for arbitrarily small tunneling. In another regime, characterized by weak interaction, but high occupation numbers, we observe a dynamical instability in the superfluid excitation spectrum. The new ground state is a superfluid, forming a 2D slab, localized along one spatial direction that is spontaneously chosen.


Introduction
Systems of ultracold atoms in optical lattices provide a unique playground for controlled realizations of many-body physics [1,2]. For sufficiently deep lattices, the kinetics is exhausted by tunneling processes, and an initially weak interparticle interaction eventually becomes important with respect to the kinetics, when the lattice is ramped up. A consequence of this competition is the quantum phase transition from a superfluid of delocalized bosons to a Mott insulator where the particles are localized at minima of the lattice by a repulsive contact interaction [3]. This effect has been observed in seminal experiments with ultracold rubidium atoms in a cubic lattice [4]. It is described quantitatively by means of the simple Bose-Hubbard model [3,5], whose parameters are the interaction energy U for each pair of particles occupying the same lattice site, and the matrix element J for tunneling between neighboring sites. Intriguing Hubbard-type physics can also be observed if the above scenario is extended to fermions, mixtures of several particle species, exotic lattice geometries, or long-ranged dipolar interaction [1,2,6].
In this paper, we consider a different type of extension of the bosonic Hubbard model, becoming relevant when the interaction between the particles is enhanced, e.g, by means of a Feshbach resonance. As long as the interaction is weak compared to the lattice potential, a system of ultracold atoms can be described, to a good approximation, in terms of the lowest-band single-particle Bloch or Wannier states, the latter being localized at the minima of the lattice [7]. Under these conditions, the Hubbard interaction U and tunneling parameter J are given by respective matrix elements with respect to the single-particle Wannier states. This approximation corresponds to degenerate perturbation theory up to first order with respect to the interaction, with only the intraband coupling induced by the interaction taken into account. However, if the interaction is stronger, higher-order corrections start playing a role. One may still describe the system in terms of lattice-site occupation numbers n j , but the occupied Wannier-like orbitals will have admixtures from higher bands, depending on the occupation. The most significant effect of the repulsive interaction will be a broadening of the Wannier-like orbitals with increasing occupation, effectively enhancing J and decreasing U. In terms of the Hubbard description, we take this into account by replacing J and U by functions Jn i ,n j and Un i of the number operatorsn i . Quantitative consequences of this kind of modification to the plain bosonic Hubbard model have been studied by several authors at a theoretical level [8,9,10]. Considering an interactioninduced modification of the Wannier functions, also additional Mott-insulator phases have been predicted [11,12]. In Ref. [13], the effect of the interaction-induced coupling to the first excited band on the Mott transition was considered. Re-entrant behaviour in the superfluid-Mott transition has also been predicted due to the interaction-induced modification of Hubbard parameters [13,14]. The effect of interaction on the tunneling dynamics in one-dimensional double-well and triple-well potentials have been studied in Refs. [15,16] where the authors found enhanced correlated pair tunneling near the fermionization limit. Moreover, occupation-number-dependent on-site interaction has been observed experimentally in the coherent dynamics of an atomic ensemble [17]. Similar occupation-dependent effects have been observed in Bose-Bose [18] and Fermi-Bose mixtures [19,20,21], and -in the latter case -have been explained theoretically in terms of occupation-dependent parameters U and J [22].
In the present work, we show that new quantum phases can arise in Hubbard models with number-dependent parameters. After writing down the effective single-band Hamiltonian including the effect of the site occupation, we find that for strong enough interaction (characterized by the s-wave scattering length a s ), there is a transition from a Mott state with one particle localized at each lattice site to a superfluid of pairs extended over neighboring sites, rather than to a superfluid of single atoms. This feature is novel, considering the fact that the extended pairs emerge in the singlespecies repulsive bosonic system without the presence of any long-range interaction. For even higher interaction strengths, the n = 1 Mott state becomes unstable towards a superfluid of pairs that are localized on single sites. Moreover, the n = 2 Mott state becomes unstable towards pair fluctuations already for very low tunneling amplitudes. Finally, we consider the regime where interaction effects are important not because of large scattering lengths, but rather because of large site-occupations numbers. In this limit, starting from the Bogoliubov approach to the homogeneous system, we find a phonon instability at a critical filling fraction. Above that fraction, the new ground state is a Bose condensate with the particle density being localized along one spatial direction that is chosen spontaneously.
This paper is arranged in the following way: In section 2, we introduce the occupation-dependent Bose-Hubbard model. In section 3, we start discussing the properties of this model. Namely, we study the instability of the Mott-insulator phase with respect to simple particle and hole excitations, leading to the usual single-particle superfluidity. In section 4, we then investigate the instability of the Mott phase with respect to the excitation of bond-centered pairs of particles being extended over neighboring lattice sites. We show that this mechanism will eventually become relevant when the s-wave scattering length is increased, and that one finds a phase transition to a superfluid of extended pairs. In section 5, proceeding to even stronger interaction, the instability of the Mott phase towards a superfluid of site-centered pairs is discussed. In this regime, moreover, the Mott insulator at a filling of two particles per site can disappear completely. Finally, in section 6, we focus on the limit where interactioninduced orbital effects play an important role because of large filling. We find that, with the increasing superfluid density, the condensate may become dynamically unstable.

The Bose-Hubbard model
The Hamiltonian in the presence of a periodic potential with lattice constant a, given by V per ( r) = V 0 [sin 2 (πx/a) + sin 2 (πy/a) + sin 2 (πz/a)], reads with bosonic field operatorsψ, mass m, and interaction strength g = 4π 2 a s /m, where a s is the s-wave scattering length. To derive a Hubbard-type description, the field operatorsψ( r) are expanded in terms of Wannier-like orbitals φ i ( r,n i ) = φ( r − R i ,n i ) localized at the lattice minima R i , namelyψ( r) = ib i φ( r − R i ;n i ) with bosonic annihilation and number operatorsb i andn i =b † ib i . Note that the "wave function" φ i depends on the number operatorn i in order to take into account interaction-induced occupation-dependent broadening. Keeping only on-site interaction, we arrive at the effective single-band Hamiltonian where and we have introduced the chemical potential µ to control the particle number. We would like to mention that in the presence of an optical lattice for high interactions the pseudo-potential form of contact interaction can still be used, when a modified scattering length which is different from the bare scattering length is applied [23,24,25,26]. In order to estimate the occupation number dependence in a mean-field way, we make a Gaussian ansatz for the Wannier-like wave functions, φ( r − R i ; n i ) = exp(−( r − R) 2 /d 2 (n i )), where the width d(n i ) is a variational parameter depending on the particle number n i , and minimize the Gross-Pitaevskii energy functional. The idea to use the width of the Wannier function as a variational parameter has also been used in Refs. [27,28,29]. Taking into account the full lattice potential (i.e., not employing a quadratic approximation for the lattice minima), for a given n i this leads to /π for the width of φ in the limit V 0 ≫ E R , where E R = π 2 2 /2ma 2 denotes the recoil energy. Note that Eq. (4) has a solution only as long as V 0 /E R ≫ d 2 (n i )/d 2 0 . Using the variational result, the tunneling parameter between two adjacent sites can be approximated by .
We would like to point out that when calculating the tunneling strength, the Gaussian approximation generally results in a lower value than the exact calculation; the exact Wannier orbital has an exponential tail which decays slower than a Gaussian. Nevertheless, our simple approximation provides us with reasonable numerical values and with a suitable model for the occupation dependence of tunneling in the regime treated here. This allows us to get a qualitative understanding of the physics at work. For the number-dependent on-site interaction strength the variational result gives The single-particle tunneling term arising from the non-on-site contributions of the quartic interaction term in Eq. (1) is exponentially smaller than J(n i , n j ) by approximately a factor of exp(−π 2 V 0 /E R /4)a s /a. Similarly the pair tunneling term is smaller than J(n i , n j ) by approximately a factor of exp(−π 2 V 0 /E R /2)a s /a. Since we are in the limit of V 0 /E R ≫ 1, these terms are neglected in Eq. (2).

Insulator to single-particle superfluid transition
Having written down a suitable model Hamiltonian describing the regime of strong interaction, we now study the transition from the Mott insulator having on average n particles per site to a superfluid of single particles/holes. For this purpose, we use a product ansatz i |Φ i for the many-body state, with the variational coherent spin-representation state [30,31], at each site i, with occupation number basis states |n i i . Here we only take into account states with one additional particle or hole, which in the Mott phase and close to the transition to the superfluid, where particle fluctuations are small, is sufficient. Accordingly, the variational mean-field energy is given by where H J = (n 2 + n)J n,n sin 2ψ/2 + (n + 1)J n+1,n sin 2 ψ + nJ n,n−1 cos 2 ψ, Minimizing the energy determines θ and ψ. While θ = 0 corresponds to an incompressible Mott-insulator state with an integer number of particles n per site (found within a finite interval of the chemical potential µ), the superfluid state is characterized by θ = 0 with order parameter b i ∼ sin 2θ. In the superfluid phase, the average particle number per site is characterized by ψ depending smoothly on the chemical potential. For ψ ≪ π/4, the transition to the superfluid occurs mainly via the creation of holes, while for ψ near π/2 particle creation is the main mechanism destroying the Mott phase. In the latter case, the Mott insulator becomes unstable when the energy cost of creating an additional particle at one site, namely U n+1 n(n + 1)/2 − µ, is overcome by the reduction in energy due to tunneling of that particle, which is on the order of z(n + 1)J n+1,n , with coordination number z = 6 for the cubic lattice. Thus, when E ss minimizes for non-zero θ, the Mott state becomes unstable with respect to single particle and hole excitations. For interaction strength a s /a = 0.15 and n = 1 , this happens at the black lines (solid or dotted) in the plane spanned by µ/V 0 and J 0,1 /V 0 in Fig. 1.

Superfluidity of extended (bond-centered) pairs
So far, we have described the usual scenario of the Mott phase becoming instable with respect to particle and hole delocalization, as it is also found for non-number-dependent Hubbard coupling J and U. However, we will now show that-as a consequence of occupation dependent hopping and on-site interaction-the Mott insulator with n = 1 can become unstable with respect to the creation of pairs of particles already before the creation of single particles becomes favorable. Consider a pair excitation with one additional particle at site i and another one at the neighboring site j, corresponding to . Such a bond-centered or extended pair excitation at ij can tunnel coherently to a neighboring bond, say ik , with k = j being another neighbor of i. Generally, bonds are considered neighbors if they share a common site. Such a pair tunneling processes occurs in second order with respect to singleparticle tunneling via the virtual site-centered pair state |P i ≡ 1 , which has larger energy. According to second-order degenerate perturbation theory, the amplitude of the pair tunneling process is given by J eff = 6J 2 2,2 /(3U 3 − 2U 2 ). On the same footing, perturbation theory gives the binding energy of −2J eff for the bondcentered pair due to number fluctuations within the pair. For a cubic lattice of sites, the bond-centered pair excitations live on an exotic lattice of coordination number z ′ = 10, being a generalization of the two-dimensional checkerboard lattice (see the rightmost drawing in Fig. 2) to three dimensions. This allows the pair to reduce its energy by 10J eff when delocalizing. In contrast, two additional particles, not forming a pair, can reduce their energy by 2 × 6 × 2J 1,2 when delocalizing on the cubic lattice of sites (coordination number 6). Thus, according to perturbation theory, the formation of a bond-centered pair is favorable if −(10+2)J eff > 24J 1,2 . For certain scattering lengths a s , this condition can be fulfilled, since the Wannier-broadening with increasing scattering lengths leads to an increase of both J 2,2 /J 1,2 and U 2 /U 3 . In such a situation, the Mottinsulator state becomes unstable with respect to the creation of bond-centered pairs rather than with respect to the creation of single-particle excitations. This happens when the delocalization energy −10J eff overcomes the energy 2(U 2 − µ) − 2J eff needed to create a pair excitation. It is interesting to note that an equivalent scenario does not happen for hole excitations, since hole excitations decrease the occupation number and with that the tunneling amplitudes. To evaluate the boundary of the n = 1 Mott-insulator phase within mean-field approximation, we construct a model for the excited bond-centered pair excitations. When the number of pairs is small compared to the number of sites, the Hamiltonian for the pairs living on top of a Mott state with one particle per site can be written as Here, L = i, j labels the bonds of the cubic lattice and LL ′ denotes pairs of nearest neighbors of these bonds as they are described by the three-dimensional checkerboard lattice (cf. Fig. 2). Moreover, we have defined the bosonic creation and destruction operators for bond-centered pair-excitationsp † L andp L , with number operatorn L = p † Lp L . As a consequence of the diluteness assumption, we have neglected the interaction between pairs, arising if pairs occupy neighboring bonds. Since the transition to a pairsuperfluid will happen with the creation of a single pair, this approximation will not influence the phase boundary. The energy of a condensate of bond-centered pairs can now be estimated in a similar fashion as before by making a product ansatz L |Φ p L of coherent states being a superposition of zero and one pair at each bond, The order parameter of the pair condensate is defined by p L = 1 2 sin(2θ e ). According to this ansatz, the variational mean-field energy per site can be written as where z ′ = 10 is the coordination number of the three-dimensional checkerboard lattice. The mean-field approach gives the same phase boundary for the appearance of a pair condensate with finite order parameter p L as the perturbation theoretical considerations of the previous paragraph. The equivalence of both approaches is generally given for an ansatz like (11) which includes only two states per site. In Figure 1 we plot the results of minimizing E ss , E ep with respect to θ, θ e for a s /a = 0.15. The stable Mott region with respect to single particle-hole excitation is given by the interior of the black solid and dotted line characterized by θ = 0. On the right hand side of the blue dashed line in Fig. 1 one finds a region where min[E ep ] < min[E ss ] with θ e = 0. Thus, here the system is characterized by p L = 0 and b i = 0, i.e., the state is a superfluid of extended pairs (ePSF).
Condensates of extended pairs have also been proposed in the context of dimer models of reduced dimensions, describing frustrated magnets like SrCu 2 (BO 3 ) 2 [32]. By approximating triplet excitations as hard-core bosons, the authors of Ref. [33] argue that for correlated hopping these bosons can condense in pairs. Such pairing processes also bear resemblance to molecular condensation due to Feshbach resonances in an optical lattice [34].
We would like to point out that triple, quadruple or higher order excitations do not play a dominant role. The effective tunneling matrix element of such excitations will be very small since it appears in third or higher order perturbation theory only. Therefore inside a n = 1 phase, triple and higher excitations cannot lower their energy efficiently by delocalization. We can, thus, exclude a superfluid of triples or higher order objects. However, there is another possible and competitive scenario we would like to mention. Instead of exciting a triple or quadruple, one can create a huge cluster of extra particles, i.e., a big spatial domain with doubly occupied sites. In this case, within each cluster, the energy of the additional particles (on top of the n = 1 Mott background) is not lowered by delocalization, but rather by the attractive interaction between them as it appears in second order perturbation theory. In the bulk of such a cluster, this gives a binding energy of −6J eff per extra particle. In comparison, in the pair superfluid each particle can lower its energy by J eff because of binding and further by another 5J eff because of delocalization (i.e., Bose condensation). Accordingly, in leading order a superfluid of bond-centered pairs on top of the n = 1 Mott insulator is equally favorable as a phase separated state with spatial domains hosting either a Mott insulator of filling n = 1 or n = 2. As a consequence, we cannot reliably exclude phase separation by means of simple variational arguments.
Before moving on, let us briefly discuss another issue: In this article, we are working in a situation with the chemical potential fixed rather than the particle number. This approach is actually quite suitable for the description of experiments with ultracold atoms, provided the atoms are trapped by a sufficiently shallow potential. In such a situation, the local density approximation applies and different regions in the trap correspond to different values of the chemical potential. However, if the trap is too steep for the local density approximation to be valid, it might introduce also new physics. Consider the following example: The phase separated state described in the preceding paragraph might not be favored in the homogeneous system. But, because it is energetically very close to the pair superfluid, it can be favored already when a slight potential difference is introduced, helping to form n = 2 Mott domains in the region of slightly lower potential energy. Such a scenario can spoil the local density approximation already for a very weak trapping potential.

Superfluidity of local (site-centered) pairs
Now, let us consider a regime that can be achieved if the interaction strength a s /a is increased further. Considering again the n = 1 Mott insulator, for increasing interaction a site-centered pair excitation, described by |P i , eventually becomes more favorable than the bond-centered excitations described by |P ij . This happens when the ratio U 3 /U 2 is reduced so much that the potential energy 3U 3 needed to create a pair of particles on the same site equals the potential energy 2U 2 required to create a pair of particles on neighboring sites. Such a situation is possible as can be derived from Eq. (4). In the limit of large V 0 ≫ E R and a s /a we can write d(n)/d 0 ≈ (gn i ) 1/5 resulting in 3U 3 − 2U 2 ≈ −0.02U 0 . If |3U 3 − 2U 2 | becomes comparable to or smaller than J 2,2 , a bond centered pair excitation |P ij can transform to a site-centered pair excitation |P i by a single-particle tunneling process described by the matrix element J pair = √ 6J 22 . In this regime, the pairs occupy the lattice given by both the sites and the bonds of the cubic lattice (see Fig. 2, center). By delocalizing on this lattice, a pair can reduce its kinetic energy by 12J pair . As long as this energy is bigger than the kinetic energy reduction 24J 1,2 which two non-paired particles can achieve by delocalization, the pair is stable towards breaking; this is the case for J 2,2 > 3/2J 1,2 . Thus, the binding mechanism of the pair is based solely on the delocalization of its center of mass. At |3U 3 − 2U 2 | ≈ 0 (given, e.g., for a s /a ≈ 0.21 when V 0 /E R ≈ 16), the n = 1 Mott insulator becomes unstable with respect to pair creation when 12J pair exceeds 3U 3 − 2µ. It is fascinating to observe the emergence of exotic lattice geometries as illustrated in Fig. 2 as a consequence of pair creation.
If the scattering length is increased further, such that 2U 2 − 3U 3 ≫ J 2,2 , sitecentered pair excitations |P i will be created rather than bond centered ones |P ij . The site-centered pair excitations can then tunnel from site to site coherently via the occupation of a virtual bond-centered pair excitation. The corresponding tunneling matrix element reads J ′ eff = 6 Moreover, the pair has a binding energy of 6J ′ eff (stemming from a small perturbative admixture of the 6 neighboring bond-centered pair states). Therefore, a site-centered pair is more favorable than two single-particle excitations if 3U 3 −12J ′ eff < 2(U 2 −12J 1,2 ). If this condition is fulfilled, the Mott insulator becomes rather unstable towards the creation of site-centered pair excitations than to the creation of single particles. The instability occurs when 12J ′ eff reaches 3U 3 − 2µ. As before, a mean-field calculation leads to the same phase boundary. We plot the boundary of the n = 1 Mott phase for a s /a = 0.3 in Fig. 3. The instability towards the creation of single particles is hardly important. It is predominantly the creation of single holes or site-centered pairs of particles which destroys the Mott phase.
Note that in the limit of U 1 ≫ J 0,1 metastable repulsively bound pairs of ultracold bosons have been observed in optical lattices [35,36]. Also, two-species mixtures of bosons with inter-species attraction trapped in an optical lattice have been shown to give rise to superfluidity of pairs [37]. In the context of dipolar atoms in a two-leg ladder, when no tunneling is present between the two legs, pair superfluidity arises due to attraction between the dipolar atoms between the two legs of the ladder [38,39]. Also, using a state-dependent optical lattice potential, it is possible to create correlated tunneling of on-site pairs, which in turn gives rise to superfluidity of local pairs [40,41]. In our present study, we find that such local pairing can emerge due to the strong occupation-dependence of tunneling and on-site interaction.
After having studied the boundaries of the Mott-insulator phase with one particle per site, let us have a look at the n = 2 Mott state. In the limit of vanishing tunneling, a Mott state with two particles localized at each site is favorable for U 2 < µ < 3U 3 − U 2 . The upper border of this interval is given by the potential energy difference of having three and two particles at a site. This difference can, in fact, become lower than the potential energy difference U 2 between two and one particle per site marking the lower border. This is the case if 3U 3 − 2U 2 < 0; then the n = 2 Mott-insulator phase is never stable with respect to the creation of particle-hole pairs, irrespective of the tunneling strength; it ceases to exist. The disappearance of the n = 2 Mott insulator coincides with site-centered pair excitations becoming more favorable than bond-centered ones in the limit of vanishing tunneling. Note that the Mott-insulator phases with higher filling, n ≥ 3, do not disappear for large interaction a s /a within the Gaussian approximation. The reason why these phases do not share the fate of the n = 2 Mott insulator is that the broadening on the Wannier-like site-wave functions φ i in response to adding one particle to that site becomes less pronounced with increasing occupation: However, one should have in mind that for strong interaction, sites occupied by three and more particles suffer strong dissipation due to three-body collisions [42,43].
One might ask about the nature of the system's ground state at fixed filling n = 2 and for 3U 3 − 2U 2 < 0, when there is no n = 2 Mott phase. At vanishing tunneling, the ground state is highly degenerate consisting of all Fock-states having occupation n i = 1 on half of the sites and occupation n i = 3 on the others. Alternatively, one might say that on top of a n = 1 Mott insulator, half of the sites are occupied by additional site-centered pairs. For small but finite hopping this degeneracy will be lifted. One can think of three possible scenarios: (i) The pairs gather in one region in space; this corresponds to a phase segregation between the n = 1 and the n = 3 Mott phases.
(ii) The pairs delocalize to form a superfluid. (iii) The pairs form a checkerboard-type insulator avoiding pairs on neighboring sites. In order to decide this question, we write down an effective Hamiltonian for the site-centered pairs: with bosonic pair annihilation and creation operatorsĉ i ,ĉ † i , and where we assume a hard-core constraint (ĉ † i ) 2 = 0. The nearest-neighbor repulsion present in the last term, with ∆ = 2 J 2 3,3 6U 4 +U 2 −6U 3 , stems from super-exchange processes between neighboring pairs. This model can be mapped to a spin-1/2 XXZ model with the first term corresponding to the XX coupling and the last one to the Z-coupling. Since (J ′ eff − ∆) ≤ J ′ eff is always true, the system will neither form the checkerboard pattern (iii) (corresponding to an antiferromagnetic state for the XXZ-magnet) nor show phase segregation (i) [44]. The system forms a superfluid of site centered pairs (ii).

Weakly interacting limit
Finally, we investigate the limit where interaction effects are important not because of a large scattering length but because of large site occupation, i.e., a s /a ≪ 1, but the mean number of particles per site n 0 ≫ 1. We assume small on-site number fluctuations δn ≪ n 0 , i.e., U n 0 /(n 0 J n 0 ) ≪ 1. In this limit, we can write the modified Hubbard Hamiltonian as where β = 3 5 π 2 and δn i,j =n i − n 0 . Here, we would like to point out the similarity of Hamiltonian (14) to the Quantum Ablowitz-Ladik (AL) model for q-deformed bosons [45], given by where [B i , B † i ] = exp[−2γN i ], and Q = 1 − exp[−2γ]. In the limit of γ → 0 and γN i ≪ 1, Eq. (18) reduces to the occupation-dependent modified Hubbard model Eq. (14) with α = γ and U n 0 = 0. It is found that in one and higher dimension the AL model contains localized solutions [46,47]. To investigate this possibility, we first solve Eq. (14) in the superfluid limit, where the order parameter reads b i = √ n 0 . To look Figure 4. Phonon velocity c as a function of the superfluid occupation number n 0 . We find that after a critical occupation number, the phonon velocity becomes imaginary, denoting a dynamical instability. The fixed parameters are a s /a = 0.01 and V 0 /E r = 10.

Conclusion and Outlook
In this paper we have predicted various effects resulting from interaction-induced band mixing in systems of ultracold bosonic atoms in optical lattice potentials. We have derived the modified bosonic Hubbard model (2) having occupation-number-dependent parameters. This model comprises an effective interaction-induced broadening of the Wannier-like single-particle orbitals, and, thus, captures also the situation when the s-wave scattering length becomes comparable to the lattice spacing, a s /a → 1. Using this model, we find that for scattering lengths a s ∼ 0.15a and lattice depths V 0 ∼ 12E R , the n = 1 Mott-insulator state can become unstable towards a superfluid which consists of bond-centered pair excitations. This scenario is novel considering the fact that the extended pairs emerge due to the occupation dependence of both the tunneling strength and the on-site interaction. For even higher interaction, the nature of the superfluid pair excitations (destroying the insulator) changes. The pairs can now occupy both the bonds on the lattice (i.e., two neighboring sites) or its sites; in that way an exotic lattice geometry as shown in the center plot of Fig. 2 emerges. Increasing the interaction further, eventually the pairs live only on the sites of the lattice. In this regime of high interaction strength, the n = 2 Mott state gets completely destroyed by the site-centered pair fluctuations. We have also looked into the regime where interaction induced Wannier-broadening arises from large filling n ≫ 1 at small scattering lengths, a s ≪ a. In this limit, we found that the superfluid becomes dynamically unstable due to the attractive nature of the occupation-dependent tunneling. The system then transforms from a uniform superfluid state to an asymmetric state which is localized in one direction and extended in the other two directions.
In future studies, we would like to study the role of dissipation in these systems. Also, a more accurate determination of the number dependence of the Hubbard parameters J n i ,n j and U n i will be required for a quantitative description of the effects described here. Finally, it would also be worth studying in detail the role of a trapping potential, as it is present in experiments.