Multi-orbital and density-induced tunneling of bosons in optical lattices

We show that multi-orbital and density-induced tunneling have a significant impact on the phase diagram of bosonic atoms in optical lattices. Off-site interactions lead to density-induced hopping, the so-called bond-charge interactions, which can be identified with an effective tunneling potential and can reach the same order of magnitude as conventional tunneling. In addition, interaction-induced higher-band processes also give rise to strongly modified tunneling, on-site and bond-charge interactions. We derive an extended occupation-dependent Hubbard model with multi-orbitally renormalized processes and compute the corresponding phase diagram. It substantially deviates from the single-band Bose–Hubbard model and predicts strong changes of the superfluid-to-Mott-insulator transition. In general, the presented beyond-Hubbard physics plays an essential role in bosonic lattice systems and has an observable influence on experiments with tunable interactions.


Introduction
Hubbard models are extremely successful in describing a variety of systems ranging from electrons in solids to bosonic quantum gases in optical lattices. While in solids all kinds of complex lattice systems are realized, optical lattices offer simple geometries in combination with precisely controllable experimental conditions. The striking idea that ultracold bosonic atoms in optical lattices constitute an excellent tool for studying the superfluid (SF)-to-Mottinsulator (MI) transition [1] has initiated an intensive investigation of bosonic lattice models. It has caused a large number of theoretical studies on the Bose-Hubbard model ranging from mean-field treatments [2][3][4] and quantum Monte Carlo [5] to density-matrix renormalization group (DMRG) [6] approaches. In addition, the influence of disorder [2,[7][8][9][10][11] and the effect of confining potentials have been discussed [12][13][14][15]. However, the first-order corrections to the Bose-Hubbard model itself, namely higher orbital tunneling (figure 1(a)) and bond-charge interactions ( figure 1(b)), have so far been mostly neglected in these approaches. Bond-charge interactions known from solids [16][17][18][19] constitute a density-induced tunneling process caused by the scattering of particles on neighboring sites. This extension to the Hubbard model was discussed in the context of superconductivity [16] and ferromagnetism [19]. While for fermionic systems these effects are partly suppressed due to the Pauli principle, they have a large impact on bosonic systems. The unique experimental access in optical lattices allows for a detailed investigation of this interaction-assisted tunneling. This density-dependent hopping was recently discussed also for boson-fermion mixtures [20,21], where the Bose-Fermi-Hubbard model does not fully cover effects arising from the interspecies interaction [20][21][22][23][24]. In this Hopping processes beyond the standard single-band Hubbard tunneling: example processes for (a) multi-orbital tunneling, (b) single-orbital bond-charge assisted hopping and (c) multi-orbital bond-charge hopping. (d) Initial and final states for one-particle tunneling from the left to the right site. The interaction causes an admixture of higher orbital states requiring a renormalization of the bare lowest-band tunneling processes.
context, the role of higher orbitals in the tunneling processes [20][21][22] was investigated. For bosonic systems, modifications of tunneling and on-site interaction due to higher bands were studied by variational mean-field methods [25][26][27][28] and by numerical exact methods restricted to either double-or triple-well systems [29][30][31][32] or to the on-site interaction [33][34][35][36]. However, in the correlated regime present in the vicinity of the MI transition, the applicability of meanfield methods to describe higher-band processes is doubtful. Thus, the challenging problem is to find a comprehensive description including both bond-charge hopping [16-20, 37, 38] and higher band processes [25-28, 30, 31, 33-36, 39, 40] that is valid also for strongly correlated systems. The recent experimental progress in accessing the occupation-dependent on-site interactions [35,[41][42][43] and the tunneling matrix element [24] necessitates an accurate theoretical method for calculating these parameters. Furthermore, the SF-MI transition in experiments with tunable interactions [23,42] is directly affected by bond-charge tunneling as well as multi-orbital extensions to conventional and bond-charge tunneling.
Here, we show that interaction-induced processes cause two substantial modifications of the Bose-Hubbard model for bosonic atoms in three-dimensional (3D) optical lattices. Firstly, the so-called bond-charge or density-induced tunneling combining interaction and hopping ( figure 1(b)) is a notable contribution to tunneling. Secondly, higher bands are occupied due to the interaction of particles on a lattice site. This increases the tunneling significantly as the tunneling in higher bands (figure 1(a)) is enhanced strongly. We present a general scheme to calculate renormalized tunneling processes by effectively incorporating higherband contributions. On this basis, we derive an effective occupation-dependent model where The latter corresponds to the background scattering length a s = 100 a 0 of 87 Rb at a lattice spacing of a = 377 nm [23]. The phase boundaries are plotted for the Bose-Hubbard model (black) and the full Hamiltonian equation (1) using meanfield (blue) and Gutzwiller (GW) (dashed red). The green line corresponds to the occupation-independent Hamiltonian equation (2). the renormalization solely depends on the solution of the on-site problem. The well-known phase diagram for bosons in lattices describing the SF-MI transition is drastically altered. The presented results have direct relevance for optical lattice experiments, which can therefore serve as an ideal testing ground for beyond Hubbard physics.
After presenting the resulting phase diagram in section 2, we discuss the extensions of the Hubbard model individually starting with the effect of bond-charge interactions in section 3. The multi-orbital renormalization of the operators for conventional tunneling, bond-charge tunneling and on-site interaction is elaborated in sections 4-7. The basic idea is to separate the strongly correlated on-site problem from the inter-site dynamics. Firstly, we solve the onsite problem numerically exactly using the Wannier states of the lattice. Subsequently, we derive renormalized matrix elements for the multi-orbitally dressed operators. Finally, we calculate the phase diagram by treating the many-site problem with mean-field theory in section 8. Note that the multi-orbital renormalization can be used for various methods to treat the many-site problem and that the results are expected to be similarly affected.

Phase diagram and model Hamiltonian
We now first explain as one central result the phase diagram for bosonic atoms (figure 2), which is derived and further discussed afterwards. The black lines depict the Mott lobes predicted by the single-band Bose-Hubbard model applying mean-field theory [4]. Our calculations (colored lines) take into account the corrections to the tunneling as mentioned above as well as multiorbital modifications of the on-site interaction. The phase diagrams are obtained for bosonic atoms of arbitrary species in cubic sinusoidal optical lattices with V (r) = V 0 i cos 2 (πr i /a), where i = {x, y, z}, a is the lattice spacing and V 0 is the lattice depth. We predict that the phase area of the MI is considerably reduced even for moderate interaction strengths and substantially deformed for stronger interactions. This corresponds to a significant shift of the critical lattice depth of the SF-MI transition (figure 3) due to an effectively increased tunneling and reduced on-site interaction. The discrepancy can reach more than 5 E R for n = 3 particles per lattice site. We expect the shift to be observable in experiments with a filling n 2 and strong or tunable interactions.
The mentioned extension to the Hubbard model leading to the phase diagram in figure 2 (red and blue lines) can be included in an effective Hamiltonian with multi-orbitally renormalized creation and annihilation operatorsb † i andb i , respectively, and the operatorñ i =b † ib i corresponding to the number of particles on site i. This Hamiltonian has the same structure as the Bose-Hubbard Hamiltonian [1,2] but with multi-orbitally renormalized tunneling matrix elements J tot n j ,ñ i and on-site interactions U MÕ n i that explicitly depend on the site occupation n i . In the following, we calculate these occupation-dependent parameters with exact diagonalization and derive the phase boundaries using a mean-field (MF) and a GW approach. Furthermore, we present a simplified Hamiltonian with occupationindependent parameters (green lines in figure 2), which is in good agreement with the full model above. This Hamiltonian is restricted to conventional and bond-charge tunneling in the lowest bandĤ with single-band operatorsb † i andb i , the Hubbard tunneling J , lowest-band bond-charge interaction J BC and renormalized on-site interactionŨ 2 for two particles, as is elaborated in 6 detail below. In general, the results for Hamiltonian (2) show that the lowest-band bondcharge tunneling is an important contribution and its influence increases with the number of particles per site. Note that the mean-field and the GW approach fully coincide for occupationindependent models (green and black lobes), whereas there is a minimal discrepancy for the Hamiltonian (1) with occupation-dependent tunneling. In the mean-field calculation, the decoupling of the lattice sites requires the use of J tot n i ,n i as an approximation for J tot n i ±1,n i and J tot n i ,n i ±1 (see appendices F and G).

Bond-charge interactions
The two-particle interaction of bosonic particles in the lowest band of the lattice is given bŷ and the Wannier functions w i (r) = w (0) (r − r i ) of the lowest band, which are maximally localized at site i. Here, a δ-interaction potential with repulsive interactions g = 4πh m a s > 0 is assumed. One of the key ingredients of the Hubbard model [1,2,44] is to restrict the twoparticle interaction to the dominating term, namely the on-site interaction U = U iiii . However, some of the neglected off-site processes correspond to scattering accompanied by hopping to a neighboring site (e.g.b † ib † jb jb j ) known as bond-charge tunneling (figure 1(c)). It is essential that this process constitutes a hopping process which physically modifies the overall tunneling in the system. Thus, although this density-dependent hopping contribution is small in comparison with the on-site interaction U , it can be neglected only if also small in comparison with the conventional tunneling J . Considering two neighboring sites i and j, the density-induced tunneling operator takes the form and J BC > 0. Due to the structure of this operator, we can combine it with the conventional tunneling J to an effective hoppinĝ As a consequence, the total tunneling increase is of the order of several tens of per cent for standard 87 Rb conditions (see figure 4, dashed line) and even more for stronger interactions and deeper lattices. Note that density-density off-site interactions (U i ji j ) [6] and correlated pair tunneling (U ii j j ) are considerably weaker than bond-charge interactions (see appendix C). The nature of bond-charge interactions can be intuitively illustrated as an effective potential. In equation (5), the termn i +n j − 1 corresponds to the density ρ BC (r) = n i |w i (r)| 2 + (n j − 1)|w j (r)| 2 on sites i and j excluding the hopping particle. By inserting explicit expressions for the integral J BC and the tunneling J , we can write the effective hopping operator (5) asĴ We can now identify the expression V (r) + gρ BC (r) as an effective tunneling potential, which is illustrated in figure 4(c). Since the density of the repulsively interacting particles is localized to the centers of the lattice sites, the effective potential corresponds to a shallower lattice, i.e. increased tunneling. Note that the effective potential can be used to determine the total tunneling via band structure calculations (see appendix E).

Multi-orbital tunneling
While the extension to the Hubbard model above is still within the lowest-band approximation, the following discussion covers the influence of higher orbitals. The calculation of the effective tunneling can be performed by taking two adjacent sites of the lattices into account (figure 1(d)). Writing initial and final states as a product of wave functions of the two sites allows us to express the renormalized tunneling in terms of single-site expectation values. Consequently, the renormalized parameters depend only on the solution of the single-site problem. The presented scheme allows us to derive the effective Hamiltonian (1) independently of the applied method for treating the on-site many-particle problem. Here, the orbital occupations are computed by exact diagonalization of n particles on a single site using a finite many-particle basis with a high-energy cut-off. The short-range interaction between the atoms is modeled as a box-shaped interaction potential (see appendix A for details). The single-particle orbitals are represented by the Wannier functions w (α) (r) with a 3D band index α = (α x , α y , α z ).
In the following, we consider multi-orbital hopping processes, where one atom tunnels from the left site (L) with initially n L particles to the right site (R) with initially n R particles ( figure 1(d)). As a direct consequence, the total tunneling J n L ,n R becomes intrinsically occupation dependent. The multi-orbital hopping operatorĴ where tunneling takes place only between equal orbitals due to orthogonality relations. It is possible to reduce this complex multi-orbital operator to an effective hopping operator The occupation-dependent tunneling matrix elements are given by J MO n L ,n R ∝ F |Ĵ MO | I , where I = (n L , n R ) denotes the initial and F = (n L − 1, n R + 1) the final state ( figure 1(d)). This expression can be separated into matrix elements for the respective lattice sites where (n) is the single-site wave function. This can be evaluated using the coefficients obtained by exact diagonalization as elaborated in appendix B. The tunneling in higher bands J α is drastically enhanced compared with the lowest band tunneling J 0 . Therefore, the contributions of higher-orbital processes are significant and grow exponentially with the lattice depth, although the population of higher bands is usually smaller than 1%. As an example, the matrix element J MO 3,3 for three bosons per site is plotted in figure 4 (green lines). In general, the effective tunneling is increased except for the elements J MO n,0 = J MO 1,n−1 where the tunneling is slightly reduced. In an intuitive picture, the many-particle calculation can be reduced to the occupation of single-particle orbitals n (α) . The fraction of particles occupying higher orbitals α tunnels via the respective tunneling matrix elements J α , as indicated in figure 4(d). In this approximation, the effective tunneling can be estimated by the weighted sum (see appendix B) This simplified approach is in qualitative agreement with the fully correlated calculation above.

Multi-orbital bond-charge interaction
It proves necessary to apply the same multi-orbital treatment of the normal tunneling developed in the last section to the bond-charge-assisted hopping. By analogy, we can define an effective multi-orbital bond-charge hoppinĝ with an occupation-dependent parameter. It depends on the single-site matrix elements (n − 1)|b

Effective tunneling
The total tunneling consists of normal and bond-charge-assisted tunneling, both effectively including higher orbital processes. The total tunneling J tot n j ,n i = J MO n j ,n i + (n i +n j −1)J BC,MO n j ,n i (12) and the individual contributions are shown in figures 4(a) and (b). The occupation-dependent enhancement of the tunneling is depicted in figure 5. The deviation from the Hubbard tunneling J can easily reach 30% for three particles per site at moderate a s and V 0 . The data in this figure are calculated using a box-shaped interaction potential (see appendix A) with a width W = 5 nm and a lattice constant a = 377 nm. It is important to note that the results are almost independent of W , provided that W is much smaller than the lattice constant. The error bars in figure 5 show the results for W = 25 nm and W → 0. Whereas the normal tunneling is increased by the multi-orbital renormalization, the density-induced tunneling is reduced. Coincidentally, both multi-orbital corrections compensate for each other and the overall deviations can be approximated surprisingly well by the singleband bond-charge tunneling for moderate particle numbers (dashed line in figure 4). This justifies the applicability of the model Hamiltonian (2), where tunneling is restricted to the lowest band, i.e. J tot n j ,n i ≈ J + (n i + n j − 1)J BC .

Occupation-number-dependent on-site interaction
In a similar fashion as the inclusion of higher orbitals alters the tunneling, also the on-site interaction U is modified. The orbital degree of freedom decreases the on-site interaction [34][35][36] as the particles tend to avoid each other. The results for U MO n , representing the eigenvalue of the single-site exact diagonalization [35], are depicted for occupation numbers n = 2-4 in figure 5. Note that for strong interactions, it may be necessary to use realistic interaction potentials. It is convenient to express U MO 2 within equation (2) by a fitted empirical functionŨ 2 /U = λ 1 + (λ 2 + λ 3 a s /a) √ V 0 /E R + λ 4 e λ 5 a s /a , which describes the on-site interaction for two particles but also serves as an approximation for three and four particles 2 .

Mean-field/Gutzwiller phase diagrams
The SF-MI phase diagram (figure 2) for the extended Hamiltonian (1) is derived within meanfield theory [2,4]. By introducing the SF order parameter ψ = b i = b † i and by neglecting the fluctuations of quadratic order, the lattice sites are decoupled (see appendix F). Within secondorder perturbation theory for the Mott state with n particles, we find an analytical expression for the SF-MI transition point withJ BC = J BC,MO n,n /J MO n,n and the unperturbed energy Due to the decoupling in the mean-field approach, only tunneling processes with n i = n j = n can be accounted for. Thus, this approach cannot fully implement the occupation-numberdependent Hamiltonian (1). This can be achieved by using the GW approach, where the energy functional is minimized with respect to the coefficients f n of the trial wave function |G = i n f n |n i (see appendix G). The GW results are depicted in figure 2 (red lines) and show nearly perfect agreement with the mean-field calculation (blue lines). The mean-field solution (13) can also be applied to the simplified occupation-independent model (2) using the lowest-band tunneling parameters J andJ BC = J BC /J as well as the fitted on-site interactionŨ 2 .
Considering the drastic simplifications of this model, the predicted Mott lobes (green lines) are in compelling agreement with the results of the occupation-dependent Hamiltonian (1).

Conclusions
We have presented an occupation-dependent model Hamiltonian incorporating multi-orbitally renormalized tunneling, density-induced tunneling and on-site interaction. The effectively increased tunneling and reduced on-site interaction cause a substantial modification of the Bose-Hubbard phase diagram for bosonic atoms. In addition, we have discussed a simplified occupation-independent model that includes the lowest-band bond-charge tunneling and reproduces the phase diagram well. In general, we have derived a multi-orbital renormalization procedure for one-and two-particle processes. It is capable of describing strongly correlated optical lattice systems accurately that cannot be calculated correctly by means of variational mean-field treatments [25][26][27][28]. Furthermore, it is important to include both bond-charge tunneling and multi-orbital corrections as extensions of the Hubbard Hamiltonian. Several effects such as the reduction of the bond-charge tunneling cannot be described with an effective 11 wave function. The presented results can be transferred to fermionic systems, quantum gas mixtures, low-dimensional systems, long-range interactions and different lattice geometries. In these systems, the occupation-dependent tunneling may also lead to novel quantum phases. Related work on defining a basis with renormalized operators can be found in [45].

Appendix A. Multi-orbital interactions
The treatment of multi-band Hamiltonians for lattices represents a challenging and still open problem posed by its complexity. A common approach is the restriction to the lowest band, which neglects interaction-induced admixtures of higher orbitals. In order to account for orbital degrees of freedom, we focus on two neighboring sites of the lattice, where we apply a multiorbital diagonalization, the so-called configuration interaction method, to each site individually (see figure 1(d), main text). In this approach, off-site interactions are neglected as the occupation of higher orbitals is mainly caused by the local on-site interaction. The 3D Wannier functions constitute the single-particle orbitals, where r L/R = (± 1 2 a, 0, 0) is the coordinate of the left and the right site, respectively, a is the lattice spacing and is the 3D band index. The maximally localized Wannier functions are obtained by means of band-structure calculations for a separable cubic lattice potential V (r) = V 0 i cos 2 (πr i /a) with the lattice depth V 0 . To include correlations of particles fully, a many-particle Fock basis |N = |n 0 , n 1 , . . . with n particles is used, where n i is the number of particles in orbital i. The Hamiltonian for a single lattice site readŝ wheren (α) =b (α) †b(α) ,b (α) † creates andb (α) annihilates a particle in the Wannier orbital α with single-particle energies (α) . In general, the ground state of the single-site problem with n particles is given as a superposition of Fock states with real coefficients c N (n). Since the Hamiltonian (A.3) preserves parity the ground state consists of many-particle states with even parity. The lowest eigenvalue of the diagonalized Hamiltonian matrix directly corresponds to the occupation-dependent on-site interaction U MO n as discussed in the main text. 12 We model the short-range interaction between the atoms as a box-shaped interaction potential with the Heaviside step function θ and a variable box width W . The influence of the finite interaction range depends only on the ratio W/a and increases with the lattice depth since the wave functions are contracted. However, for W a the results (see figure 5, main text) are almost independent of W . The interaction integrals in equation (A.3) are given by whereg is fixed by the condition U (0000) (g, W ) = g d 3 r |w(r)| 4 with g = 4πh m a s and s-wave scattering length a s . Note that the δ-interaction potential can be understood as the limiting case of the box interaction potential. However, it behaves differently if including an infinite number of orbitals requiring a regularization procedure as described in [33] (for numerical diagonalization see [46]). For real interaction potentials, the mathematical subtlety of the δ-interaction is not present and under certain experimental conditions [35] decay processes in higher orbitals may even effectively lead to a finite Hilbert space. Here, for the numerical diagonalization, nine bands are accounted for in each spatial direction restricted to the energetically lowest 6 × 10 3 × n 2 many-particle states.

Appendix B. Multi-orbital tunneling
In the following, we consider hopping processes, where one atom tunnels from the left site (L) with initially n L particles to the right site (R) with initially n R particles (see figure B.1(d), main text). The tunneling operator accounting for all possible orbital hopping processes readŝ L with tunneling matrix elements Because of orthogonality relations, only matrix element with α = β have nonzero contributions. Thus, the tunneling between the left and the right site simplifies tô with matrix elements J α ≡ J αα , which only depend on the band index α x in the x-direction. The respective tunneling matrix elements J α for the band α x are plotted in figure B.1, which are positive for even bands and negative for odd bands. The tunneling amplitude is strongly enhanced for higher orbitals due to the considerably larger overlap of the Wannier functions. As the occupation of higher orbitals is decoupled from off-site interactions and tunneling in our model, it is possible to reduce the complex multi-band hopping to an effective hopping J MO n L ,n R . It describes the transition from an initial state | I to the final state | F with As the initial and final states (B.3) are product states, the expression separates into two terms for the left-and the right-hand sites Using the many-particle ground state (A.4) obtained by exact diagonalization, we can define Note that since the coefficients c N can be chosen to be real, the conjugated term reads (n − 1)|b (α) † | (n) = j * (α) n−1 = j (α) n−1 . The effective tunneling can thus be written as Due to parity conservation of the Hamiltonian (A.3), the ground states at the left and the right site are superpositions of Fock states |N (n) with even parity. Tunneling in an odd orbital would alter the parity of both sites, resulting in final states with vanishing coefficients. Consequently, tunneling can only occur in even orbitals α with positive tunneling matrix elements J α . It turns out that the effective tunneling is increased by the tunneling in higher bands with the exception of the processes J MO n,0 = J MO 1,n−1 . Here, the lack of interaction in the left or the right site prohibits tunneling via higher bands. The interaction at one site, however, depopulates the Fock state |n, 0, 0, . . . , which results in slightly lowered tunneling energy. The process J MO 1,0 , where no interactions are present, is equivalent to the uncorrelated tunneling J = J 0 .
Intuitively, the mean-field wave functions of the left and right sites can be used to calculate the multi-orbital tunneling [28]. The coefficients can be obtained from the orbital occupation numbers c (α) n = n (α) n / √ n, with the expectation value n (α) n = (n)|n (α) | (n) . This allows us to estimate the multi-orbital by (B.11)

Appendix C. Off-site interactions
The full lowest-band interaction with interaction integrals (3) is commonly restricted to the on-site interaction. However, off-site interactions between neighboring sites can reach the same order of magnitude as the tunneling. Expanding the lowest-band Hamiltonian accordingly leads to three distinct physical processes, namely bond-charge interaction (density-dependent hopping), correlated pair tunneling and density-density interaction [16][17][18][19]37]. While the latter process has to be compared with the onsite interaction, the other two constitute tunneling processes. The full lowest-band Hamiltonian with nearest-neighbor interaction readŝ with matrix elements U = U iiii , J BC = −U iii j , J pair = U ii j j /2 and V = U i ji j . In figure C.1 these parameters are plotted as a function of the lattice depth. The offsite density-density interaction is very small compared to the on-site interaction and can consequently be neglected. This also applies for the correlated pair tunneling, which is even negligible when compared with the single-particle tunneling matrix element J . The bond-charge tunneling matrix element, however, reaches 10% of the conventional tunneling amplitude for intermediate and deep lattices. In addition, it scales with the total particle number and can thus be a significant contribution to the total tunneling. Note that the multi-orbital renormalization can influence the individual parameters very strongly.

Appendix D. Multi-orbital bond-charge tunneling
In analogy to the multi-orbital tunneling, the orbital degree of freedom also affects the bondcharge interaction. The single-band bond-charge or density-induced tunneling is introduced in the main text (see equation (3) and figure 1(b)). As discussed before, we can restrict the calculation to two neighboring sites and write the multi-orbital bond-charge operator (see figure 1(c), main text) aŝ where the Wannier functions are chosen to be real and the sign P = (−1) α x +β x +γ x +δ x is parity dependent. Note that for a δ-interaction potential, we have J BC αβγ δ = J BC αβδγ . The effective multiorbital bond-charge tunneling readŝ Note that (n − 1)|b (β) †b(γ ) †b(δ) | (n) = j (δγβ) n−1 and that because of the same arguments as in the last section, the orbitals α and β + γ + δ need to be even; thus the sign P in equation (D.1) vanishes.
the effect of the bond-charge interaction as a normal tunneling within the effective potential V eff = V (r) + g(ρ(r) − |w j | 2 ) experienced by the hopping particle (see figure 4(c), main text). However, since U iii j = U i j j j , one may also defineρ i j = (n i − 1 2 )|w i | 2 + (n j − 1 2 )|w j | 2 . In this case, the effective potential for a homogeneous filling n 1 can be written as V (r) + g i (n − 1 2 )|w i | 2 . This form of the potential can be used to perform a band structure calculation in the effective potential [21]. Note, however, that in general the effective potential is not separable into its spatial directions. Furthermore, in a band structure calculus the Wannier functions are adapted to the effective potential whereas in equation (E.5) the Wannier functions of the bare lattice potential, are used. The perturbation series up to third order in ψ reads with the second-order correction E 2 (n) = n|J BC n,n (n − 1) + 1| 2 E 0 (n) − E 0 (n − 1) + (n + 1)|J BC n,n n + 1| 2 E 0 (n) − E 0 (n + 1) + 1.
The boundary of the SF-MI phase transition is given by the Landau criterion E 2 (n) = 0 for second-order phase transitions, which can be solved for µ (see figure 1, main text).
with z = 6 nearest neighbors. The system is in an MI state when all fluctuations vanish, i.e. when |G = i |n i is the energetically lowest state. The resulting GW phase diagram is plotted together with the mean-field results in figure 1 in the main text. Both the approaches, GW and mean-field, are in almost perfect agreement.