Twonniers: Interaction-induced effects on Bose-Hubbard parameters
Mark Kremer Rashi Sachdeva Albert Benseny Thomas Busch
cond-mat.quant-gas

	

	


Institut für Physik, Universität Rostock, Albert-Einstein-Str. 23,18059 Rostock, Germany

rashi.sachdeva@oist.jp
Quantum Systems Unit, Okinawa Institute of Science and Technology Graduate University, Okinawa 904-0495, Japan

Quantum Systems Unit, Okinawa Institute of Science and Technology Graduate University, Okinawa 904-0495, Japan

Quantum Systems Unit, Okinawa Institute of Science and Technology Graduate University, Okinawa 904-0495, Japan




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.


67.85.Hj, 67.85.-d

Twonniers: 
 Interaction-induced effects on Bose–Hubbard parameters
    Thomas Busch
    December 30, 2023
===================================================================






§ 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 <cit.>. 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 <cit.>. 

Theoretically, ultracold atoms in optical lattices can be described by a Bose–Hubbard model <cit.>, 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 <cit.>.

However, while mathematically convenient, single particle Wannier functions neglect certain physical effects, such as the broadening of the localized wave functions 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 <cit.> or when using optical lattices for metrology <cit.>, 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 <cit.>. 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  <cit.>.  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. <ref> we provide a brief review of the conventional way of calculating the Hubbard parameters using the single-particle Wannier function approach.
Then, in Sec. <ref> 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. <ref> to calculate the parameters of the modified Bose–Hubbard Hamiltonian,
which are interpreted in Sec. <ref> in comparison to those obtained from single-particle Wannier functions.
Finally, we discuss possible applications and conclude in Sec. <ref>. 



§ THE BOSE–HUBBARD MODEL
 


The starting point for our analysis is the Hamiltonian for a Bose gas, given by 

    Ĥ = Ĥ_SP + Ĥ_I ,
 
where the single-particle term includes the kinetic energy and the optical lattice potential,

    Ĥ_SP =	∫d Ψ̂^† () [-ħ^2/2m∇^2+V_L()]Ψ̂ ().

Here m is the atomic mass.
The term including the point-like interactions is given by

    Ĥ_I = g/2∫d Ψ̂^† ()Ψ̂^† ()Ψ̂ ()Ψ̂ () ,

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(), and bosonic  annihilation and creation operators, â_i and â_i^†, for each lattice site as

    Ψ̂ ()=∑_i f_i()â_i
    	  with  ∫d  	f_i^*()f_j()=δ_ij.

A convenient and common choice for the orthonormal functions in a lattice potential are the well-known Wannier functions <cit.>, which are localized at the individual lattice sites. 
The single-particle Wannier function at lattice site i in the Bloch band α is defined as 

    w^α_i() = w_i,x^α(x) w_i,y^α(y) w_i,z^α(z) ,

and the components in each direction can be written in terms of the Bloch functions ϕ^α_k(x) as 

    w_i,x^α(x) = 1/√(N_x)∑_ke^-i k x_iϕ^α_k(x),

where N_x is the number of lattice sites along the x-direction (equivalent expressions exist for the other spatial directions), and x_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 single-particle 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

    J=∫d  w_i^*() [-ħ^2/2m∇^2+V()]w_i(),

where only the nearest-neighbor overlaps are taken into account, and the interaction part of the Hamiltonian leads to the onsite interaction amplitude

    U=g∫d  |w_i()|^4 .





§ 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. (<ref>)).
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()f_i() by two-particle Wannier functions, but leave terms of the form f_i()f_j() (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(), interacting via a point-like potential. The Hamiltonian is given by

    Ĥ = ∑_k=1^2[-ħ^2/2m∇^2_k+V_L(_k)]+g/2δ(_1-_2),

and its corresponding delocalized eigenfunctions Φ_j(_1,_2)  can be used as a basis to construct the localized  (two-particle) functions

    W_i(_1,_2)=∑_j c_jΦ_j(_1,_2)   with ∑_j|c_j|^2=1.

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 <cit.>

    M_i=∫d_1d_2 W_i^*(_1,_2)(_1^2+_2^2)W_i(_1,_2).

This allows us to define the single-particle single-site densities from the two-particle wave functions as
|W_i(,)|. In order to fulfill the orthogonality condition in Eq. (<ref>) this density needs to be normalized as

    ∫d |W_i(,)|!=1,

which also assures the fulfilment of the particle statistics,

    [a_i,a_j^†] = δ_i,j  and  [a_i,a_j]=[a_i^†,a_j^†]=0.



To compare the single particle and two-particle Wannier functions, we show in Fig. <ref> their respective densities computed in a one-dimensional  potential V_L(x)=V_0sin^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.









§ 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. (<ref>), takes the form

    Ĥ_I   =g/2∑_ijkl∫d  f_i^*()f_j^*()f_k()f_l()â_i^†â_j^†â_k â_l 
       =1/2∑_ijklU_ijklâ_i^†â_j^†â_k â_l.

As we are only interested in the ground state, the Wannier functions and the two-particle wave functions based on Eq. (<ref>) 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

    f_i()f_j()
    
    |W_i(,)|    if i=j, 
    
    w_i()w_j()    if i≠ j,

which should be compared to the standard way of calculating using single-particle Wannier functions

    f_i()f_j()
    
    w_i()w_j()
      ∀ i,j.

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 (<ref>), 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

    U_iiii = g∫d  f^4_i() ,

which under the substitutions of Eqs. (<ref>) and (<ref>) becomes

    U_iiii^W   =g∫d  |W_i(,)|^2 ,
    
    
    	U_iiii^w   =g∫d |w_i()|^4 .


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

    U_iijj = g∫d  f^2_i()f^2_j() ,

and become after substitution

    U_iijj^W   =g∫d  |W_i(,)|  |W_j(,)| ,
    
    
    	U_iijj^w   =g∫d |w_i()|^2 |w_j()|^2 .


The next effect is associated with terms including â_i^†â_j^†â_i â_j, and it can be interpreted as two indistinguishable processes:
the interaction between particles 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

    U_iiij = g∫d  f^3_i()f_j() ,

which, after the substitutions, become

    U_iiij^W   =g∫d  |W_i(,)| |w_i()| |w_j()| ,
    
    
    	U_iiij^w   =g∫d |w_i()|^3 |w_j()| .
 




§ 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
g_1D=-2ħ^2/ma_1D,
with a_1D=-d_⊥^2/2a_s(1-Ca_s/d_⊥),
where C ≃ 1.4603 and d_⊥=√(2ħ/mω_⊥)  <cit.>. 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. <ref>. 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. <ref>(c)), where 
U_iiii^W is smaller than U_iiii^w. 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 potential 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. <ref>(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. <ref>(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. <ref>(e)). 


These results are consistent with the situation where the interaction strength is changed while keeping the lattice depth constant (see Fig. <ref>).
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. <ref>(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. <ref>(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 <cit.>, which take the two-particle tunnelling and the cross tunnelling terms into account <cit.>.





§ POSSIBLE APPLICATIONS AND CONCLUSIONS





To summarize, we have calculated the parameters for the Bose–Hubbard model by consistently including on-site 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. <cit.>, the atomic clock shift of ^87Rb 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 <cit.>. When calculated using single particle Wannier functions, the onsite interaction term is independent 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. 





Greiner:02 M. Greiner, O. Mandel, T. Esslinger, T.W. Hänsch, and I. Bloch, Nature 415, 39 (2002).


Yukalov:09 V. Yukalov, Laser Phys. 19, 1 (2009).


lewenstein_book M. Lewenstein, A. Sanpera, and V. Ahufinger, Ultracold Atoms in Optical Lattices, (Oxford Univ. Press, 2012).


Jordens:08 R. Jördens, N. Strohmaier, K. Günter, H. Moritz, and T. Esslinger, Nature 455, 204 (2008).


Fisher:89 M.P.A. Fisher, P.B. Weichman, G. Grinsten, and D.S. Fisher, Phys. Rev. B 40, 546 (1989).


Jaksch:98 D. Jaksch, C. Bruder, J.I. Cirac, C.W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).


vanOosten:01 D. van Oosten, P. van der Straten, and H.T.C. Stoof, Phys. Rev. A 63, 053601 (2001).


Zwerger:03 W. Zwerger, J. Opt. B: Quantum Semiclass. Opt. 5, S9 (2003).


Campbell:06 G.K. Campbell, J. Mun, M. Boyd, P. Medley, A.E. Leanhardt, L.G. Marcassa, D.E. Pritchard, and W. Ketterle, Science 313, 649 (2006).


latticeclock:14 B.J. Bloom,	T.L. Nicholson, J.R. Williams, S.L. Campbell, M. Bishof, X. Zhang, W. Zhang, S.L. Bromley, and J. Ye, Nature 506, 71 (2014).


Li:06 J.B. Li, Y. Yu, A. Dudarev, and Q. Niu, New. J. Phys. 8, 154 (2006).


Liang:09 Z.X. Liang, B. Hu, and B. Wu, arXiv:0903.4058 [cond-mat.str-el] (2009).


Schneider:09 P.-I. Schneider, S. Grishkevich, and A. Saenz, Phys. Rev. A 80, 013404 (2009).


Luhmann:12 D.S. Lühmann, O. Jürgensen, and K. Sengstock, New. J. Phys. 14, 033021 (2012).


Zhu:15 S. Zhu and B. Wu, Phys. Rev. A 92, 063637 (2015).

wan37 G.H. Wannier, Phys. Rev. 52, 191 (1937).
kohn73 W. Kohn, Phys. Rev. B 7, 4388 (1973).


Stollenwerk:11 T. Stollenwerk, D.N. Chigrin, and J. Kroha, J. Opt. Soc. Am. B 28, 1951 (2011).
olshanii98 M. Olshanii, Phys. Rev. Lett. 81, 938 (1998).



rashi1 R. Sachdeva, S. Johri, and S. Ghosh, Phys. Rev. A 82, 063617 (2010).

rashi2 R. Sachdeva and S. Ghosh, Phys. Rev. A 85, 013642 (2012). 


benseny A. Benseny, J. Gillet, and Th. Busch, Phys. Rev. A 93, 033629 (2016).


reshodko I. Reshodko, A. Benseny, and Th. Busch, Phys. Rev. A (2017); arXiv:1703.02189 (2017).



ebhm:lewenstein M. Maik, P. Hauke, O. Dutta, M. Lewenstein, and J. Zakrzewski, New J. Phys. 15, 113041 (2013).


Cornell:02 D. M. Harber, H. J. Lewandowski, J. M. McGuirk, and E. A. Cornell, Phys. Rev. A 66, 053616 (2002).
Baym:96 G. Baym and C. J. Pethick, Phys. Rev. Lett. 76, 6 (1996).




