Interaction broadening of Wannier functions and Mott transitions in atomic BEC

Superfluid to Mott-insulator transitions in atomic BEC in optical lattices are investigated for the case of number of atoms per site larger than one. To account for mean field repulsion between the atoms in each well, we construct an orthogonal set of Wannier functions. The resulting hopping amplitude and on-site interaction may be substantially different from those calculated with single-atom Wannier functions. As illustrations of the approach we consider lattices of various dimensionality and different mean occupations. We find that in three-dimensional optical lattices the correction to the critical lattice depth is significant to be measured experimentally even for small number of atoms. Finally, we discuss validity of the single band model.


I. INTRODUCTION
Numerous many-body phenomena have been recently demonstrated with Bose-Einstein condensates (BEC) in optical lattices [1,2,3]. Number squeezing has been observed with 87 Rb atoms in a one-dimensional lattice of pancake-shaped wells [1], and superfluid to Mottinsulator transitions have been witnessed with such atoms in three-dimensional and one-dimensional optical lattices [2] . Such transitions were predicted by theoretical studies based on the Bose-Hubbard model [4] and by microscopic calculations of the model parameters for BEC in optical lattices [5,6].
Very important question is whether it is possible to observe superfluid to Mott-insulator transitions for the mean occupation number n larger or even much larger than one? Phenomenological single band Bose-Hubbard model indeed predicts such transitions. Previous calculations of the model parameters J, hopping amplitude, and U , on-site interaction, were based on the lowest band Wannier functions for a single atom in an optical lattice. Repulsive interaction between the atoms for n > 1 may cause the wave function in each well to expand in all directions, not only affecting the on-site interaction U [7] but also strongly enhancing tunneling J between neighboring wells. This is especially significant in lower dimensional lattices with transverse potential bigger than the lattice wells where large occupations can be achieved without substantial three-body collisional loss. In order to provide theoretical guidance for experimental observation of Mott transitions in such systems, it is very important to obtain accurate critical parameters of the lattice potential for lattice occupations beyond unity.
Here we show how to construct an orthogonal basis of Wannier functions with mean-field atomic interactions taken into account. We use it to obtain renormalized values of parameters J and U , from which critical depth of the potential V c is calculated for various lattices of different dimensionality and mean occupation. For the cubic optical lattice with n = 2 or larger, our result is no-ticeably larger than that calculated without taking into account interaction. This increase is more pronounced for the anisotropic cases with stronger lattice potentials in one or two directions. For the case of one-dimensional lattice of pancake-shaped wells [1] or two-dimensional lattice of tubes [3], our results are several times larger than critical values calculated from one-atom Wannier functions. This is in agreement with the experimental findings that much higher lattice potentials are needed to reach the transition point in such cases.
Kohn developed variational approach to calculate electronic Wannier functions in crystals [8]. We modify this procedure by minimizing on-site energy self-consistently taking into account interaction between atoms.
In the last section we address validity of the singleband Bose-Hubbard model constructed with variational Wannier functions. The conditions for the model to be valid need to be modified from those for a single particle case since the interaction between the particles alters the band structure substantially [9]. For the model to be valid two conditions have to be fulfilled: (i) when the number of particles in a well changes by one the variational Wannier function should not change significantly and (ii) collective excitations of the atoms within each well should be less energetically favorable than atom hopping between the wells.

II. BOSE-HUBBARD MODEL AND WANNIER FUNCTIONS
For bosonic atoms located in the lattice potential V (r) and described by boson field operators ψ(r), the Hamiltonian field operator is ∇ 2 + V (r) ψ(r) where a s is the atoms' scattering length and m is the mass. To illustrate our methods we use as an example isotropic cubic lattice. We assume that the boson field operator may be expanded as ψ(r) = i b i W (r − r i ), where b i is the annihilation operator for an atom in the Wannier state of site r i . Substituting this expansion into the Hamiltonian we obtain a problem of lattice bosons. We consider the case when the number of atoms per cite n i fluctuates around average number n. This results in the standard Bose-Hubbard Hamiltonian where the effective on-site repulsion U , the hopping amplitude J and the on-site single-atom energy I are defined by where g = 4πa s 2 /m and a is the lattice vector. On-site energy f is defined as with the bare on-site interaction U 0 We assumed that the Wannier function does not change much for small fluctuations of the number of atoms. Offsite interactions are also neglected. In case of more than one atom per site the presence of other atoms does modify the Wannier function of an atom. Below we describe our strategy for its selfconsistent calculation. We start with a trial wave function localized in each well, g(r − r i ). A Wannier function corresponding to the lowest Bloch band may be constructed according to Kohn's transformation: where the integral is over the first Brillouin zone and For an odd Wannier function, the cosine function should be replaced by the sine function. One can show that such Wannier functions are normalized and are orthogonal to each other for different wells. We vary the trial function to minimize the on-site energy f [10]. We note that another method to calculate the Wannier functions including interaction effects self-consistently may be used for small interactions. Starting with nonlinear time-independent Gross-Pitaevskii equation one may calculate periodic Bloch states u k (r) defined as by expanding them in Fourier series and solving nonlinear system of equations. Then, a set of Wannier wave functions for the band in question is defined by This procedure fails for large interactions because the bands develop loops and become not single-valued [9].

III. SUPERFLUID TO MOTT-INSULATOR TRANSITIONS
We consider three optical lattice systems which are relevant to experiments: (i) isotropic three-dimensional optical lattice, (ii) anisotropic three-dimensional lattices, and (iii) the situation when the lattice potential is present only in one or two directions and confinement in other directions is provided by relatively weak harmonic trap. Following standard practice, we will use the lattice period π/k, atomic mass m, and recoil energy E r = 2 k 2 /2m as the basic units.
Three pairs of counter-propagating laser beams with wavelength 2π/k propagating along three perpendicular directions create potential Isotropic cubic lattice is created by the beam of equal intensity. In this case V Anisotropic cubic lattices can be created by choosing intensity of one or two beam to be much large than other.
Below we study the case when ω ⊥ ≫ µ, where µ is the chemical potential of the atoms, thus the weak optical lattice is effectively one-dimensional or two-dimensional and transverse motion is frozen to the ground state of the transverse confinement.
Transverse motion can also be decoupled in the experimentally relevant case when the lattice potential is present only in one or two directions and atoms are confined in other directions by relatively weak harmonic trap: V T (r ⊥ ) = 1 2 mω 2 ⊥ r 2 ⊥ . According to existing experiments, in our calculations through this work, we choose the 87 Rb atoms in F = 2, m = 2 state with scattering length a s = 5.8 nm and the laser wavelength of 852 nm for the three-and twodimensional lattices and 840 nm for the one-dimensional lattice. All numerical results are obtained using 21 lattice wells in each direction with periodic boundary condition. Convergence has been checked using 41 wells for some of the key results.
In each case we calculate parameters of the Bose-Hubbard model based on the variational approach described in the previous section. The critical condition for superfluid to Mott-insulator transition has been found approximately as where z is the number of the nearest neighbor sites [11]. By substituting the parameters into the critical condition, we can map out the critical potential strength as a function of mean occupation.
In the following, we report our findings for isotropic and anisotropic three-dimensional lattices, onedimensional lattice of pancake wells, and two-dimensional lattice of tubes.

A. Isotropic cubic lattice
In the case of isotropic cubic lattice we choose variational trial function to be in the form g(r) = g(x)g(y)g(z), with g(u) = (1 + αu 2 )e −u 2 /σ 2 , where α and σ are variational parameters. Then the Wannier function must also be of the product form W (r) = w(x)w(y)w(z), with the one-dimensional functions w(u) and g(u) related by the one-dimensional version of Kohn's transformation. All the three-dimensional integrals in Eq. (2)-(15) can then be reduced to one-dimensional ones, greatly simplifying the calculations.
Our calculations proceed as following. For a given V 0 and n, we start with certain initial parameters α and σ to obtain a trial Wannier function through Kohn's transformation and calculate the on-site energy f . The procedure is repeated by varying the parameters until the on-site energy f is minimized. The resulting variational Wannier function will depend on both n and V 0 . If only the on-site single-atom energy I is minimized, one obtains the single-atom Wannier function W 0 (r) which only depends on V 0 . We find that interaction broadens Wannier functions, as a result U 0s is always larger than U 0 , but we also notice that effective interaction U can be larger  (7) respectively. The derivative in (3) is calculated by Chebyshev fitting to function f . Interaction parameter U0S calculated with single particle Wannier function is defined as U0S = g dr|W0(r)| 4 , where W0(r) is a single-atom Wannier function. than U 0 (see Fig. 1). So phase transition is more complex than we expected.
Once the Wannier function is determined, we can calculate the Bose-Hubbard parameters U and J. In Fig. 3, we depict the ratio U/zJ (z = 6) as a function of the mean occupation n for several values of the potential strength V 0 . The decreasing trend can be understood as following. The total interaction energy increases with n, making the Wannier function broader, hence the interaction parameter U becomes smaller, J proportional to overlap between neighboring Wannier functions becomes larger, and as a result the ratio decreases. The intersection with the line of critical condition (in Fig. 3 the line with positive slope obtained from Eq. (15)) then yields the mean occupation for which these potentials are crit-   [12], and are also consistent with experimentally determined range for the critical potential [? ]. For n > 1, the mean field repulsion makes the critical potential noticeably higher. Starting from n = 3 the correction to the critical depth of the lattice has to be clearly observable experimentally and effects of interaction has to be taken into consideration.

B. Anisotropic cubic lattices
Our procedure can also be applied to the case of an anisotropic lattice. We model the system as a lower dimensional problem with the reduced interaction parameter g d obtained by multiplying g by the integral of |ψ ⊥ | 4 , where ψ ⊥ is the single-atom ground state wave function in a well of the transverse potential [5]. In the harmonic approximation, the wave function can be found exactly, and the reduced interaction parameter is given by V ⊥ for the quasi-one-dimensional lattice and g 2 = g π √ V ⊥ for the quasi-two-dimensional lattice. In the calculations discussed below, we take V ⊥ = 80E r .
To find the Wannier functions for the lower dimensional lattices, we use these reduced interaction parameters in our procedure, replacing all the threedimensional integrals in Eqs. (5) and (5) by lower dimensional ones. The critical lattice potential V c calculated using such variational Wannier functions is depicted in Fig. 4 for the one-and two-dimensional models. For comparison, we also include results calculated using the one-atom Wannier function. The increase of critical potential due to mean-field repulsion on the Wannier functions is somewhat bigger in the lower dimensional cases.

C. Lattices in one or two directions
BECs in one-dimensional lattice of pancake-shaped wells and two-dimensional lattice of tube-shaped wells have been studied in experiments [1? ]. Because of the large transverse dimensions of such wells, many atoms can be held in a well without suffering too much threeatom collisional loss, opening the possibility of studying superfluid/Mott-insulator transition for relatively large n [7,13]. In a theoretical investigation, Oosten et al [7] considered the interaction effect by using a transverse wave function in the Thomas-Fermi approximation without modifying the single-atom Wannier function in the lattice direction(s). Here we extend their work by considering the interaction effect on the Wannier functions as well.
For the pancake like BEC array, the transverse wave functions are approximated by the Thomas-Fermi wave function φ T F (r ⊥ ) of the BEC within the pancake plane, which is defined by for µ > V T (r ⊥ ) = 1 2 mω 2 ⊥ r 2 ⊥ and vanishes otherwise. According to the experimental data, we take ω ⊥ = 19 × 2πs −1 . We begin by writing the Wannier function in the form, W (r) = w(r L )φ(r ⊥ ), where φ is the wave function for the transverse direction(s), and w is the Wannier function in the lattice direction(s), both to be determined variationally by minimizing the on-site energy. The part of the on-site energy involving φ is just the n-particle Gross-Pitaevskii energy in the transverse potential and with the interaction parameter g modified into g d by multiplying the integral of |w(r L )| 4 . In the Thomas-Fermi approximation, this 'transverse energy' is given by f ⊥ = 2n−1 3 nmω 2 ⊥ g 1 /π for the one-dimensional case and f ⊥ = 5n−2 10 (9mω 2 ⊥ n 2 g 2 2 ) 1/3 for the two-dimensional case. The total on-site energy is the sum of this 'transverse energy' and n times of the single-atom energy of the lattice Wannier function: Lattice Wannier function, obtained by the procedure of Kohn's transformation and minimization of the on-site FIG. 6: Energy associated with hopping (process 1) has to be smaller than energy to excite the many-body state in well (process 2). Many-body excitation is schematically depicted as a single atom excitation. energy, will be affected by the interaction because the 'transverse energy' depends on it through the reduced interaction parameter g d . After w(r L ) is determined variationally, the Bose-Hubbard parameters J and U can be calculated immediately. In Fig. 7, we show the critical potential V c for the case of one-dimensional lattice with transverse trap frequency ω ⊥ /2π = 19 s −1 and 120 s −1 .
For comparison, we also show the corresponding results obtained using the single-atom Wannier function of the lattice and the Thomas-Fermi transverse wave function. It is clear that V c is raised dramatically due to the broadening of the Wannier function. In the experiment of Ref. [1], the magnetic trap potential is 19 s −1 . The transverse trap frequency is enhanced to 120 s −1 if the optical confining potential with V 0 = 50E r is turned on, and the mean occupation number is n ∼ 50. Evidence from Bragg interference pattern shows that the critical value of the lattice potential should be somewhat larger than 44E r . This observation is contradictory to the prediction based on the single-atom Wannier function, but is consistent with our result based on the variational Wannier function.
In the case of two-dimensional lattice, our results for the critical lattice potential are shown in Fig.  7(b) for ω ⊥ /2π = 24 s −1 which is used in [? ]. We predict V c ∼ 33E r for n ∼ 100, while the single-atom Wannier function yields V c ∼ 27E r . The largest lattice potential used in the experiment was 12 E r , so further experiment is needed to verify the theoretical predictions.

IV. VALIDITY OF THE SINGLE-BAND MODEL
In this section, we discuss the conditions for the single band Bose Hubbard model to be valid. First, we make general remarks and then give quantitative examples relevant for the case of the isotropic cubic lattice.
Assumption that the boson field operator may be expanded as ψ(r) = i b i W (r − r i ) requires that the Wannier functions do not change substantially when the number of atoms in a well changes by one. A good criteria for this condition to be fulfilled is that interaction energy calculated with the Wannier function does not change much when number of particles changes by one Note that the value of U can still be quite different from the one calculated with a single particle Wannier function. When the condition is fulfilled, the second condition is that the excitations within the ansatz have to be the least energetical. That is the hopping of the atoms from well has to be more energetically favorable than excitation of atoms in each well to the many-body excited state (see Fig. 6). If we consider two neighboring wells the energy of the ground state is E 0 = 2nI + U 0 n(n − 1).
The energy associated with hopping is ∆E 1 = 2nI + U 0 (n − 2)(n − 1) + n(n + 1) 2 It has to be much smaller than the energy of the first excited many-body state that we denote ∆ We plot the criteria from Eq. (18) for isotropic lattices on Fig. 7. It is much smaller than unity. To estimate the effect of many-body excitation within a single well, we neglect hopping amplitude J, since close to Mottinsulator transition it is much smaller than atom's interaction. Also for the experimentally relevant region of the potential depths the potential can be well approximated by a harmonic potential. In the harmonic potential the lowest many-body excited mode is associated with the center of mass motion -Kohn mode [15]. As a result ∆ ∼ ω. Since we neglect the tunneling we may start directly with variational form for the Wannier function in a well. We take W (x, y, z) = W (x)W (y)W (z), where W (u) = C(1 + βu 2 )e −γu 2 . Similar to previous section for a fixed V 0 and n we minimize on-site energy f . From the results shown in Fig. 8 it is clear that the single-band model is applicable in this case: the energy associated with atom's hopping is much smaller than the energy required to excite the atoms inside of the wells.