Twonniers: Interaction-induced effects on Bose-Hubbard parameters

We study the effects of the repulsive on-site interactions on the broadening of the localized Wannier functions used for calculating the parameters to describe ultracold atoms in optical lattices. For this, we replace the common single-particle Wannier functions, which do not contain any information about the interactions, by two-particle Wannier functions ("Twonniers") obtained from an exact solution which takes the interactions into account. We then use these interaction-dependent basis functions to calculate the Bose--Hubbard model parameters, showing that they are substantially different both at low and high lattice depths, from the ones calculated using single-particle Wannier functions. Our results suggest that density effects are not negligible for many parameter ranges and need to be taken into account in metrology experiments.


I. INTRODUCTION
Ultracold atoms in optical lattices have been a recent topic of significant interest, as they can be used to perform quantum simulations of fundamental models of many-body physics, which are often difficult to access using traditional condensed matter systems [1][2][3]. The perfect periodicity of optical lattices allows to mimic the crystalline environments electrons experience in solids and unprecedented control over the kinetic properties of the atoms is possible by tuning the lattice depths. Furthermore, the interaction properties between the ultracold atoms can be changed using techniques like Feshbach resonances. This has opened up many new avenues of research, particularly in the field of condensed matter and atomic physics, and made it possible to study quantum phases and quantum phase transitions over a wide range of parameters [1][2][3][4].
Theoretically, ultracold atoms in optical lattices can be described by a Bose-Hubbard model [5][6][7][8], which stems from a mapping of the continuous system to the lattice by using site localized single-particle Wannier functions. The static and dynamics properties of the gas are then described by two main parameters: the hopping term, which accounts for bosons tunneling between neighboring sites, and the on-site interaction term, which accounts for the repulsive energy when two particles sit at the same lattice site. The competition between these parameters (commonly determined by calculating overlap integrals using single-particle Wannier functions) characterizes the Mott-insulator/superfluid transition [1].
However, while mathematically convenient, single particle Wannier functions neglect certain physical effects, such as the broadening of the localized wave functions * Electronic address: rashi.sachdeva@oist.jp due to repulsive on-site interactions when two or more bosons occupy the same lattice site. This can have significant effects when trying to make precision measurements [9] or when using optical lattices for metrology [10], as the energy scales that govern the behavior of the atoms are typically small.
Recently, a number of theoretical efforts have been made to incorporate the effects of interaction on the Wannier functions using mean-field and numerical approaches [11][12][13][14][15]. In addition, there has been strong experimental evidence of the broadening of Wannier function at high fillings, when high-resolution spectroscopy showed non-uniform frequency shifts for different occupation numbers per site [9]. It is therefore important to include the effects of modified densities due to the repulsive interactions when calculating the Bose-Hubbard parameters. In this work we suggest to do this by using the exact two-particle wave functions ("Twonniers"), obtained after solving the two-particle Schrödinger equation with contact interaction. For comparison, we also perform calculations using the single particle Wannier functions. To the best of our knowledge, this is the first study where the expansion is directly performed in terms of the two-particle wave functions, which has an implicit dependence on repulsive atom-atom interactions.
Our presentation is organized as follows. In Sec. II we provide a brief review of the conventional way of calculating the Hubbard parameters using the single-particle Wannier function approach. Then, in Sec. III we introduce the two-particle wave functions that include the interaction effects by solving the two-particle Schrödinger equation with contact interaction. These wave functions are used in Sec. IV to calculate the parameters of the modified Bose-Hubbard Hamiltonian, which are interpreted in Sec. V in comparison to those obtained from single-particle Wannier functions. Finally, we discuss possible applications and conclude in Sec. VI.

II. THE BOSE-HUBBARD MODEL
The starting point for our analysis is the Hamiltonian for a Bose gas, given bŷ where the single-particle term includes the kinetic energy and the optical lattice potential, Here m is the atomic mass. The term including the pointlike interactions is given bŷ where g = 4π 2 a s /m is the interaction strength related to the s-wave scattering length, a s . The bosonic field operators,Ψ andΨ † , can be expanded into a series of orthonormal functions, f i (r), and bosonic annihilation and creation operators,â i andâ † i , for each lattice site aŝ A convenient and common choice for the orthonormal functions in a lattice potential are the well-known Wannier functions [16,17], which are localized at the individual lattice sites. The single-particle Wannier function at lattice site i in the Bloch band α is defined as and the components in each direction can be written in terms of the Bloch functions φ α k (x) as where N x is the number of lattice sites along the xdirection (equivalent expressions exist for the other spatial directions), and x 0 i is the center of the i-th trap. It is important to note that the Wannier functions are not eigenfunctions of the system and that, as singleparticle functions, they do not contain any information about possible scattering effects due to multi-particle occupancy of a site. Also, for small interaction energies the particles can be considered to be confined in the lowest Wannier orbitals because the energy separation between the lowest and first excited band is quite large compared to interaction energy. We work in this regime and from now onwards will drop the band index α.
The hopping amplitude in the Bose-Hubbard model can then be calculated as where only the nearest-neighbor overlaps are taken into account, and the interaction part of the Hamiltonian leads to the onsite interaction amplitude

III. TWO-PARTICLE WAVE FUNCTIONS
The effect of the repulsive scattering interaction depends on both the interaction strength g and the density distribution of the wave function (see Eq. (3)). Therefore, it is important to choose the correct form for the orthonormal functions with which one performs the expansion: since the interactions are local and the functions are localized the density distribution should take the interaction into account if two (or more) particles are at the same lattice site. We will therefore in the following replace terms of the form f i (r)f i (r) by two-particle Wannier functions, but leave terms of the form f i (r)f j (r) (i = j) to be described by single-particle Wannier functions.
To find the two-particle Wannier functions we solve the Schrödinger equation for two particles in a sinusoidal potential, V L (r), interacting via a point-like potential. The Hamiltonian is given bŷ and its corresponding delocalized eigenfunctions Φ j (r 1 , r 2 ) can be used as a basis to construct the localized (two-particle) functions Since the interactions raise the energies, we use the eigenfunctions of the two lowest bands. To determine the coefficients c j , we assume that the particles are well localized at each lattice site, using as the criteria for localization the minimization of the second moment [18] This allows us to define the single-particle single-site densities from the two-particle wave functions as |W i (r, r)|.
In order to fulfill the orthogonality condition in Eq. (4) this density needs to be normalized as which also assures the fulfilment of the particle statistics, The red (solid) line corresponds to the twoparticle single site density obtained after numerically solving the Schrödinger equation with the Hamiltonian (9) for 9 traps, with lattice depth V0 = 1.5Er and scattering length as = 100a0, where a0 is the Bohr radius. The blue (dashed) line corresponds to the square of single-particle Wannier functions for the same lattice parameters. The lattice depth is given in units of the recoil energy Er = π 2 2 /2ma 2 , where a is the lattice spacing of the sinusoidal optical lattice potential. The inset shows a zoom-in on the tails of the densities, clearly showing the broadening of two-particle density compared to the density of the single-particle Wannier function.
To compare the single particle and two-particle Wannier functions, we show in Fig. 1 their respective densities computed in a one-dimensional potential V L (x) = V 0 sin 2 (πx/a). One can clearly see that, as expected, the repulsive interaction leads to a broadening of two-particle Wannier function, which eventually results in significant change in the Bose-Hubbard parameters. However, one can also see that the wings of the two particle Wannier function at the position of the neighbouring lattice sites are suppressed, which is due to the orthogonality requirement between two of the modified Wannier functions.
In the next section, we use this two-particle wave function and density to construct the different terms in the Hamiltonian and compare them to the ones using only single-particle Wannier function solutions.

IV. MODIFIED BOSE-HUBBARD HAMILTONIAN
The effects of the interactions between the particles are fully contained in the interaction termĤ I , which, after inserting the expansion of Eq. (4), takes the form As we are only interested in the ground state, the Wannier functions and the two-particle wave functions based on Eq. (9) can be chosen to be real and we will therefore neglect the complex conjugates below. The parameters U ijkl can then be calculated using the substitution which should be compared to the standard way of calculating using single-particle Wannier functions Here we have introduced the labels W and w which will be used below to distinguish, respectively, terms calculated from the two-particle Wannier function density or from single-particle Wannier functions. The hopping term in the Bose-Hubbard model depends only on the single-particle Wannier functions as it comes from the non-interacting part of the Hamiltonian (2), and it is therefore is not affected by these substitutions.
To explicitly identify the different physical processes that are summarized in the interaction term, we will in the following group the different terms into four categories. The first one is the one where two particles are at the same site and interact with each other. The associated terms includeâ † iâ † iâ iâi and their corresponding amplitude is given by which under the substitutions of Eqs. (15) and (16) becomes The second group corresponds to terms with operatorsâ † iâ † iâ jâj , (i = j), which describe the joint tunneling of two particles between two neighbouring lattice sites, i.e. the particles hop together from one lattice site to another. The coupling amplitudes associated with this process are given by and become after substitution U w iijj = g dr |w i (r)| 2 |w j (r)| 2 .
The next effect is associated with terms includinĝ a † iâ † jâ iâj , and it can be interpreted as two indistinguishable processes: the interaction between particles (a) as = 100a0 . The insets shows the behavior for shallow lattices. For the numerical calculation 9 traps have been taken into account. Bottom row: Ratios of (c) Uiiii, (d) Uiijj, and (e) Uiiij, calculated with the two methods (single-particle and two-particle Wannier functions) for as = 100a0 (solid blue) and as = 400a0 (dashed red).
at neighbouring sites or cross tunneling of particles. As these processes only involve a single particle at each site, one gets U W ijij = U w ijij = U w iijj . Finally, the last effect is associated with terms includingâ † iâ † iâ iâj , which describes single-particle tunneling between an empty and an already occupied neighbouring trap. The coupling amplitudes for this process are given by which, after the substitutions, become

V. RESULTS AND DISCUSSIONS
In the following we will numerically compute and compare the interaction parameters for the single-particle and the two-particle Wannier function approach. To avoid complications from the regularized delta function in three dimensions, all calculations are done in one dimension, assuming a tight harmonic confinement of the atoms in the transverse direction (of frequency ω ⊥ ). However, all calculations are conceptually straightforward to extend to higher dimensions. Adjusting the coupling constant g to one dimension can be done via 19]. In the following we choose ω ⊥ = 2π × 10 4 Hz. The results for two different values of the scattering length (a s = 100 a 0 and a s = 400 a 0 ) and as a function of the lattice depth are shown in Fig. 2. It can be seen that the overlap integrals U iiii , which describe the on-site interaction, are generally in good agreement with each other for both approaches. The biggest deviations appear for shallow lattices (see Fig. 2(c)), where U W iiii is smaller than U w iiii . The difference stems from the fact that the repulsive interaction leads to a broadening of two particle density and consequently a reduction in its maximal amplitude, which directly translates into a smaller magnitude of the interaction coefficient for the two-particle Wannier approach. For deeper lattices, i.e. larger poten-   tial energies, the broadening is reduced and the two quantities have similar values. The crossing between U iiii and J, which is visible in the inset of Fig. 2(a), corresponds to the parameter range where tunneling starts to dominate over the interaction effects. Since at the crossing point the two relevant values of U iiii differ by about 10%, an effect on the Mott-transition point can be expected.
Similar differences between the two methods can also be noted for the overlap integrals for the correlated pair tunneling, U iijj , where for shallow lattices the integral based on the two-particle Wannier functions is larger than the one based on the single-particle functions. Here the extended size of the localised functions due to the repulsive interactions leads directly to a larger overlap between neighboring sites. On the other hand, for deeper lattices, the pair-tunneling coupling calculated from the two-particle functions becomes an order of magnitude smaller than that from the single-particle functions. This is due to the fact that even at higher lattice depths the single particle Wannier function density and the two particle density have different behaviour in their tails, although their bulk density becomes almost identical. In this regime, the magnitude of the tail of the single particle Wannier density is higher than the one of the two particle density, leading to a larger overlap between neighboring densities, and thus to higher values of U iijj (see also Fig. 2(d)). Finally, the density dependent couplings U iiij show a difference for shallow lattices, which can be explained in the same way as for the interaction terms above (see Fig. 2(e)).
These results are consistent with the situation where the interaction strength is changed while keeping the lattice depth constant (see Fig. 3). The on-site interaction and interaction-mediated tunneling terms, U iiii and U iiij , do not show much difference between the two methods, but the two-particle tunneling coupling, U iijj is much more severely affected. For a comparatively deep lattice (V 0 = 20E r , Fig. 3(b)) the two-particle tunneling amplitude calculated using the two-particle Wannier approach increases faster than the one based on the single-particle Wannier functions, and the two methods do not coincide anywhere in the plotted parameter regime. However, for a shallower lattice (V 0 = 10E r , Fig. 3(a)) a crossing can be seen, as the two curves associated to U iijj are closer together. This leads to the conclusion that the effects of the interactions can have significant influence on the parameters of the Bose-Hubbard model, and should be taken into account in particular in metrology experiments. It also provides justification for the use of extended Bose-Hubbard models [20,21], which take the two-particle tunnelling and the cross tunnelling terms into account [22][23][24].

VI. POSSIBLE APPLICATIONS AND CONCLUSIONS
To summarize, we have calculated the parameters for the Bose-Hubbard model by consistently including onsite density effects. This was done by replacing the commonly used single-particle Wannier functions by two particle Wannier functions, which result in a broadening of the density due to repulsive interactions. Given the experimental control parameter of the optical lattice depth and the scattering lengths, we have shown that in certain regimes the Bose-Hubbard parameters show substantial deviation from the results using single-particle Wannier functions and that terms such as the correlated pair tunnelling can be become important, even though they are usually neglected.
These results are hence of principle interest for current and future experiments in the field of ultracold atoms in optical lattices, especially to account for non-uniform shifts in atomic clock frequencies due to the collision of atoms. In a recent experiment by Campbell et al. [9], the atomic clock shift of 87 Rb was measured, and found to decrease with increasing number of atoms per site. Other works have also shown that the clock frequency shift is directly proportional to the onsite interaction strength [25,26]. When calculated using single particle Wannier functions, the onsite interaction term is in-dependent of the occupancy of lattice sites, and hence cannot explain the decrease of the clock shift with increasing occupancy. However, the presented technique takes into account the effect of repulsive interactions implicitly, and the resulting broadening of the two-particle single-site density and the decrease of the magnitude of onsite interaction term U iiii , can explain the decrease of clock shift.