Bose-Einstein condensation in an optical lattice: A perturbation approach

# Bose-Einstein condensation in an optical lattice: A perturbation approach

C. Trallero-Giner Faculty of Physics, Havana University, 10400 Havana, Cuba    Victor Lopez-Richard Departamento de Fisica, Universidade Federal de São Carlos, 13.565-905, São Carlos, São Paulo, Brazil    Ming-Chiang Chung Institute of Physics, Academia Sinica, Taipei 11529, Taiwan    Andreas Buchleitner Physikalisches Institut der Albert-Ludwigs-Universität, Hermann-Herder-Str. 3, D-79104 Freiburg
July 12, 2019
###### Abstract

We derive closed analytical expressions for the order parameter and for the chemical potential of a Bose-Einstein Condensate loaded into a harmonically confined, one dimensional optical lattice, for sufficiently weak, repulsive or attractive interaction, and not too strong laser intensities. Our results are compared with exact numerical calculations in order to map out the range of validity of the perturbative analytical approach. We identify parameter values where the optical lattice compensates the interaction-induced nonlinearity, such that the condensate ground state coincides with a simple, single particle harmonic oscillator wave function.

###### pacs:
03.75.Be, 03.75.Lm, 05.45.Yv, 05.45.–a

## I Introduction

The dynamics of Bose Einstein Condensates (BEC) in harmonically confined optical lattices is one of the key objects of cold matter research shin (). While the fundamental equation which describes the mean field solution, the Gross-Pitaevski equation (GPE) Groos1 (), is readily amenable to a numerical solution, explicit analytical expressions for the condensate ground state and the associated chemical potential are scarce. Variational solutions have been proposed in perez (), and quasiclassical approximations were implemented, providing several analytical results for the GPE Konotop (). Eigenenergies and eigenmodes of the non-linear Schrödinger equation in a parabolically confined optical lattice can be characterized through asymptotic expansions of Mathieu functions, in the tight-binding approximation ana ().

We will demonstrate in the present contribution that closed analytic results can indeed be derived, and that they permit not only qualitative but also quantitative insight into the competition between the interaction-induced nonlinearity on one hand, and the redistribution of the particles over the spatially modulated potential, on the other. Comparison with exact numerical results will finally allow us to demarcate the range of validity of our perturbation approach, with respect to parameter values employed in state of the art experiments.

## Ii Analytical expressions

We consider the one-dimensional nonlinear GPE in a stationary optical lattice Groos1 (),

 {−ℏ22md2dx2+12mω2x2+λ1D|Φ|2 −sERcos2(2πdx)}Φ=μΦ, (1)

where is the frequency of the confining harmonic trap Exp1 (), the atomic mass, the recoil energy Kramer (), an adjustable parameter defined by the laser intensity, the laser wavelength, and the atom-atom interaction strength. With this one dimensional ansatz, we implicitly assume that the dynamics in the transverse direction be frozen – e.g., by sufficiently strong confinement or some other adiabaticity condition Berg (). Furthermore, we imply the number of atoms per site to be large enough to justify the mean field limit Pitaevskii (); Pitaevskii2 (). The normalization condition for reads

 1=∫∞−∞dx|Φ|2. (2)

Equation (1) admits the following integral representation for trallero ():

 Φ(x)=∫∞−∞G(x,x′)f(x′)dx′, (3)

with the Green function of the linear operator , and

 f=(μ+sERcos2(2πx/d)−λ1D|Φ|2)Φ

the inhomogeneity of the differential equation .

Recall that for Eq. (3) to hold, we have to impose that all involved functions tend to zero as . Thus, on the basis of the general theory of Fredholm integral equations Mihling (); petrovskii (); trallero2 (), the harmonic oscillator functions represent a complete basis set which span the Hilbert space of the integral equation (3). Using the spectral representation of , we can rewrite the order parameter as

 Φ=∞∑n=0φn(xlo)Cn(μ), (4)

where is the harmonic oscillator wavelength. Inserting (4) into (3), the vector coefficient is restricted to obey the relation Mikhlin ()

 [Δ(μ)+ΛCτ⋅T⋅C−VoP]C=0, (5)

with the transpose of the vector . and are fourth and second rank tensors, respectively, and defined in Appendices A and B. Furthermore, , , and

 Δnm=(n+12−μℏω)δnm. (6)

### ii.1 Chemical potential

We tackle the problem under the assumption that the atom-atom interaction, the nonlinear term , and the optical lattice potential can be considered as perturbations with respect to the trap potential . In this case, the dimensionless chemical potential and the vector in Eq. (5) can be expanded in series,

 C = C(0)+C(1)+C(2)+…, (7) μℏω = μ(0)+μ(1)+μ(2)+…, (8)

where the quantities and are understood to be of the same order in and . Expansion of Eq. (5) to second order in and leads, after lengthy but straightforward manipulation, the following expression for the system’s ground state energy:

 μℏω = 12+ΛT0000−VoP00−3Λ2∞∑m=1|T0002m|22m+ (9) 4ΛVo∞∑m=1T0002mP02m2m−V20∞∑m=1|P02m|22m.

Using the properties of and given in the Appendices A and B we obtain

 μℏω = 12+Λ√2π−Vo2[1+exp(−α2)]+ (10) ΛVo√2πexp(−α2)∞∑m=1α2mm2mm!− 3Λ22π∞∑m=1(2m−1)!24m(m!)2− V2o4exp(−2α2)∞∑m=122mα4m2m(2m)!,

with the ratio of trap to laser wavelength. Finally, the series in (10) can be summed up with the help of Eqs. (27-30) in Appendix C, and we obtain for the chemical potential, at second order in and :

 μℏω = Λ√2π+1−Vo2−Vo2exp(−α2) (11) +ΛVo√2πexp(−α2){Ei(α2)−C−lnα2} −0.033106×Λ2 −V2o4exp(−2α2){Chi(2α)−C−ln2α}.

### ii.2 Order parameter

Following the same procedure as above, and thus (through (4)) the normalized order parameter can be expressed, at first order in and , as

 Φ=φ0(xlo)−∞∑m=1ΛT0002m−VoP02m2mφ2m(xlo). (12)

With the expressions derived in the Appendices A and B for the matrix elements and , this can be rewritten as

 Φ = φ0(xlo)+∞∑m=1(−1)m+12m{Λ√(2m)!√2π22mm!− (13) Vo2m−1√(2m)!α2mexp(−α2)}φ2m(xlo).

Since, according to Eqs. (31-34) of Appendix C,

 F(z,α)=exp(−α2)√π1/2∞∑m=1(−1)m2m−12m√(2m)!α2mφ2m(z)

equals

 F(z,α) = 1√loπ1/2exp(−α2)exp(−z2/2)× ∫α0[exp(y2)cos(2yz)−1]dy2y,

and, moreover, with (LABEL:f-38),

 G(z)=∞∑m=1(−1)m+1√(2m)!2m22m√2πm!φ2m(z)

can be condensed into

 G(z)=exp(−z2/2)√2π√loπ1/2∫√2/21exp(−z2y2(1−y2))−y1−y2dy,

we obtain the following closed expression for the order parameter:

 Φ=φ0(xlo)+Λ√loG(xlo)+Vo√loF(xlo,α). (14)

### ii.3 Numerical Solution

For a numerical solution of Eq. (1), we subdivide the dimensionless coordinate , , such that , , with a grid size . For , a three-point approximation with uniform spacing allows to rewrite the second derivative and the differential equation (1) as difference equations:

 (z2i2+Λ∣∣¯¯¯¯¯ϕi∣∣2−Vocos22πlozid+1δ2)¯¯¯¯¯ϕi −¯¯¯¯¯¯¯¯¯¯ϕi−12δ2−¯¯¯¯¯¯¯¯¯¯ϕi+12δ2=μℏω¯¯¯¯¯ϕi, (15)

where is the local, normalized order parameter which abides by the boundary conditions .

Equation (15) defines a tri-diagonal eigenvalue problem with the normalization condition or, equivalently, , where . Since the resulting matrix (15) involves the non-linear term , an iterative procedure is needed for a self-consistent solution. To avoid numerical instabilities trallero2 (); Pu (), we inserted the iterative wave function

 F(k)i=√(1−ε)∣∣~ϕ(k−1)i∣∣2+ε∣∣~ϕ(k)i∣∣2 (16)

in (15), at each new iteration , with , and starting condition . This procedure was repeated until and , where and are the desired accuracies for the order parameter and the chemical potential, respectively. Our iterative solutions then converge to a relative uncertainty of for and the normalized order parameter.

## Iii Results

Figure 1 displays the variation of the chemical potential, , as a function of the dimensionless atom-atom interaction . The analytical solution (10) is represented by solid, dashed, and dash-dotted lines, while symbols represent the numerical solutions of Eq. (1). The inset illustrates the confining potential , as a function of for and . In our calculation of , we fixed and checked the variation of the perturbative solution as compared to the numerical one, for several values and of the normalized wavelength. For the relative error between the numerical solution and Eq. (11) increases as increases, while for the worst comparison is obtained at . In the case of repulsive interaction, numerical and perturbative calculation exhibit small variations with . In general, Eq. (11) describes the values of the chemical potential better for repulsive than for attractive interaction. This can be understood with the following argument: In trallero (); trallero2 () it was shown shown that, at , perturbation theory matches the numerical solution very well in the interval . When the optical lattice is switched on (), we can argue that an effective renormalization of the non-linear parameter takes place. Since Eq. (11) includes a negative quadratic term in , we expect that the effective range of where perturbation theory is applicable will shift to , i.e., more repulsive or less attractive atom-atom interactions. Naturally, the behavior of the chemical potential with limits this simple description.

In Fig. 2, we plot the dependence of on , for several values of the strength of the periodic potential . In the calculation we took . The figure shows perfect agreement of Eq. (11) with the numerical result for ranging between and , while, at , we observe large discrepancies.

Figure 3 displays the normalized order parameter, for four values of , and , and fixed , . In general, we obtain a remarkably good agreement between the analytical prediction (14) and the numerical result, for all considered (attractive) values of the dimensionless interaction parameter , and witness the interplay between the optical lattice and harmonic oscillator potential: The harmonic oscillator ground state is modulated by an oscillatory behavior, induced by the optical lattice potential. We further note that the wave function gets more localized, and the maximum of increases, as the non-linear potential becomes more attractive. Panel (a) of Fig. 4 shows the variation of with the strength of the periodic potential , for the repulsive case . As increases from 0.5 to 10, the lattice-induced density oscillations manifest with increasing amplitude. A comparison between attractive and repulsive interactions, for , is given by panel (b) of Fig. 4: Since in Figs. 3 and 4, the wave function exhibits oscillations at (), which are quenched because of the asymptotic decay of the oscillator wave functions. These oscillations are illustrated in Fig. 4(b), where the solid line represents the potential (see inset of Fig. 1), along with the wavefunction for the same values of and . The positions of the observed maxima (minima) in the order parameter correspond to the minima (maxima) of the resulting effective potential .

The extrema of exhibit different behavior depending on the non-linear interaction . As the condensate becomes effectively less confined for increasing repulsive interaction, the amplitudes of the maxima and minima of are enhanced. From the physical point of view it is clear that a strongly attractive interaction (in the figure: ) should flatten the lattice-induced oscillations of the condensate. In panel (b) we compare the evolution of the order parameter from attractive to repulsive interaction. We note that, as increases, the condensate spreads, i.e., the wave function delocalizes, while its maximum decreases. The analytical solution (14) is less accurate at and (with an error less than ), in agreement with the results of Fig. 1 and the anticipated range of validity of perturbation theory. Nevertheless, in this particular case the observed agreement with the exact solution is quite remarkable at , and shows that the closed analytical expressions derived here allow a fair representation of the BEC wave function in an optical lattice. Almost no differences can be observed on the scale of the figure between numerical and analytical solutions for .

### iii.1 Quenching of the non-linear interaction

Since the repulsive atom-atom interaction and the optical lattice potential have opposite signs, the set of parameters (, , , ) can be chosen such as to minimize the effect of the non-linearity. To study this quenching effect, in the validity range of our perturbation approach, we consider the order parameter given by Eq. (14). In the inset of the left panel of Fig. 5, we monitor the dependence of and on . While the function is universal, shows a strong dependence on the parameter . Moreover, for the contribution of to the order parameter, induced by the optical lattice, resembles the function , though with opposite sign. If we impose the condition

 ΛG(x/lo)+VoF(x/lo,α)→0, (17)

then can be solely described by the harmonic oscillator wave function . The range of values of , , and that satisfy (17) is represented by the hatched region in Fig. 5. At , Eq. (17) defines the curve shown in the figure. We have noted that at the functions and reach their minimum and maximum values, respectively. In the right panel of the figure we compare the harmonic oscillator wave function with the order parameter as predicted by (14). We chose three different points , , and in the parameter space of the left panel for this comparison: Indeed, does not fulfill the condition (17), and, consequently, does not match . The case belongs to region , but not to the curve . Hence, we observe a small discrepancy between both functions, mainly at . Case , which sits right on top of the curve , yields a perfect match between harmonic oscillator wave function and order parameter. Exact numerical solution of the GPE (1) with the parameters defined by corroborates this result, as indicated by the open circles in the figure.

### iii.2 Validity of the perturbation approach

For a vanishing lattice, , our closed analytical solution for the chemical potential and its comparison with the numerical solution provides universal criteria for its range of validity: in the interval , we derive an accuracy better than trallero (); trallero2 () for Eq. (11). For nontrivial values of the lattice strength, two more independent parameters, and come into game. A necessary condition for the validity of the perturbative results requires that the functions , , and , that contribute to Eq. (11), simultaneously be smaller than one. This criterion allows to confine the range of validity of our analytical results in a 2D map spanned by and , as displayed in Figs. 6(a) and (b). There the contribution to (11) is plotted for three different values of , , and . Fig. 6(a) thus defines different vs. regions where Eq. (11) can be trusted: In region I, the quadratic term in is less important than the linear one, while, in region II, we need to make one more distinction, depending on the values. For example, if , we have that represents the main contribution to , and vice versa if , where . A similar analysis applies in regions III and IV, with respect to the relative contributions of and .

Finally, let us note that experimental data as reported in Ott () on Rb condensates were obtained for and . This working point is indicated by a solid circle in Fig. 6(b), and shows that our perturbative method is applicable in the experimental parameter range.

## Iv Conclusions

We provided closed analytical expressions for the chemical potential and the particle density of the GPE ground state in an optical lattice. These solutions were obtained by a perturbative expansion in the lattice strength and in the (attractive or repulsive) atom-atom interaction, with a well-controlled range of validity in the associated parameter space.

Interestingly, the interaction-induced non-linearity may be quenched by the presence of the lattice. Under certain conditions (see left panel of Fig. 5), we predict that the particle density of the repulsively interacting system is given by .

The solutions derived here can serve as a useful tool to study a weakly interacting BEC in a not too deep 1D optical lattice. On this basis, it is possible to develop a unified and comprehensive picture of Bogoliubov equations, the time-dependent GPE and of collective excitations (a subject under investigation). The model here developed can be generalized to two and three dimensions, and also for two component BECs.

###### Acknowledgements.
C.T-G. is grateful to the Alexander von Humboldt Foundation for financial support, and for the hospitality enjoyed during his stay at the Max-Planck-Institut für Physik komplexer Systeme in Dresden. V. L.-R. acknowledges financial support by the Brazilian agencies FAPESP and CNPq.

## Appendix A Tensor T

Using the harmonic oscillator wave function

 φn(x)=(1π1/22nn!l0)1/2exp(−x22l02)Hn(xl0), (18)

with the Hermite polynomials, the matrix elements are defined by

 Tplmn = 1π√2n+m+l+pn!m!l!p!× ∫∞−∞exp(−2z2)Hn(z)Hm(z)Hl(z)Hp(z)dz.

It is possible to perform the above integral, and we get trallero2 (); lord ()

 Tp,l,n,m = (−1)M−m−p2M−12π√2n+m+l+pn!m!l!p!× Γ(M−l+12)Γ(M−n+12)Γ(M−n−l+12)× 3F2(−m,% \ −p,−M+n+l+12;−M+l+12,−M+n+12;1),

with the gamma function, the generalized hypergeometric function Gradshteyn80 (), and . Consequently, we arrive at the following useful relations: i) , if is an odd number, ii), for Gradshteyn80 (),

 Tpln0=2s−1π2Γ(s−l)Γ(s−p)Γ(s−n)√2n+l+pl!p!n!, (21)

and, iii), for ,

 T0002m=(−1)m√(2m)!√2π22mm!. (22)

## Appendix B Tensor P

The contribution of the optical lattice is represented by the two dimensional matrix elements given by

 Pkp = 1√π2k+pk!p!∫∞−∞[cos2(2πlodz)× (23) Hk(z)Hp(z)exp(−z2)]dz.

For this integral is equal to Gradshteyn80 ()

 Pkp = 12m+1√k!(k+2m)!{δm,0+ (24) (−1)mb2mexp(−b2/4)L2mk(b22)},

where and are the Laguerre polynomials. The symmetry of the Hermite polynomials imposes that if . Using (24) the following relations hold: i) For ,

 P00=12+12exp[−(2πlod)2], (25)

and, ii), for and ,

 P02m=(−2)m2√(2m)!(2πlod)2mexp(−(2πlod)2). (26)

## Appendix C Series

We can sum up the series

 F1=32π∞∑m=1(2m−1)!24m(m!)2 (27)

noting that

 ∞∑m=1(2m−1)!23m(m!)2(x2+1)m=−ln(12√2x2+12(x2+1)+12)2. (28)

Hence, .

Furthermore, we have the series Gradshteyn80 ()

 (29)
 ∞∑m=122m2m(2m)!α4m=Chi(2α)−C−ln2α, (30)

where is the exponential integral, is the cosine hyperbolic integral, and is Euler’s constant.

For the summation of the first part of the series in (13), we note that

 F2(z,α) = ∞∑m=1(−2)mα2m2√(2m)!φ2m(z) (31) = exp(−z2/2)2√l0π1/2(expα2cos2αz−1).

The series

 G1(z,α)=12∞∑m=1(−2)mα2m2m√(2m)!φ2m(z) (32)

is related to the function through the differential equation

 dG1(z,α)dα=F2(z,α)α. (33)

The solution of (33) is given by

 G1(z,α) = 12√l0π1/2exp(−z22) (34) ∫α01y[exp(y2)cos(2yz)−1]dy.

To perform the sum of the first series in (13), we use that

 F3(z,c) = 1√l0π1/2[cexp(−c2z2+z22)−exp(−z22)] = ∞∑m=1(−1)m√(2m)!2mm!(1−1c2)mφ2m(z).

Now the series

 g(z,c)=∞∑m=1(−1)m+1√(2m)!m2mm!(1−1c2)mφ2mz) (36)

is related to the function through the equation

 dgdβ=2β(1−β2)F3, (37)

with The differential equation (37) admits the solution

 g(z,c)=2exp(−z22)√l0π1/2∫1/c1exp(−z2y2(1−y2))−y1−y2dy. (38)

## References

• (1) Th. Anker, M. Albiez, R. Gati, S. Hunsmann, B. Eiermann, A. Trombettoni, and M. K. Oberthaler, Phys. Rev. Lett. 94, 020403 (2005); Y. Shin, G.-B. Jo, M. Saba, T. A. Pasquini, W. Ketterle, and D. E. Pritchard, Phys. Rev. Lett. 95, 170402 (2005); A. V. Ponomarev, J. Madroñero, A. R. Kolovsky, and A. Buchleitner, Phys. Rev. Lett. 96, 050404 (2006); K. Winkler, G. Thalhammer, F. Lang, R. Grimm, J. H. Denschlag, A. J. Daley, A. Kantian, H. P. Büchler, and P. Zoller, Nature 441, 853 (2006); T. Roscilde and J. I. Cirac, Phys. Rev. Lett. 98, 190402 (2007); J. Brand and A.R. Kolovsky, Eur. Phys. J. D 41, 331 (2007).
• (2) E. P. Gross, Nuovo Cimento 20, 454 (1961); L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 40, 646 (1961) [Sov. Phys. JETP 13, 451 (1961)].
• (3) V. M. Pérez-García, H. Michinel, J. I. Cirac, M. Lewenstein, and P. Zoller, Phys. Rev. Lett. 77, 5320 (1996); ibid, Phys. Rev. A 56, 1424 (1997); V. I. Yukalov, E. P. Yukalova, and V. S. Bagnato, Laser Phys. 12, 1325 (2002); ibid, Phys. Rev. A 66, 025602 (2002);
• (4) V. Konotop and P. Kevrekidis, Phys. Rev. Lett. 91, 230402 (2003); T. Hyouguchi, R. Seto, M. Ueda, and S. Adachi, Ann. of Phys. 312, 177 (2004); D. Witthaut, H. J. Korsch, J. Phys. A: Math. Gen. 39, 14687 (2006)..
• (5) A. M. Rey, G. Pupillo, Ch. W. Clark, and C. J. Williams, Phys. Rew. A 72, 033616 (2005).
• (6) M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, E. A Cornell, Science 269, 198 (1995).
• (7) M. Krämer, L. Pitaevskii, and S. Stringari, Phys. Rev. Lett. 88, 180404 (2002).
• (8) K. Berg-Sørensen and K. Mølmer, Phys. Rev. A 58, 1480 (1998); S. Burger, F. S. Cataliotti, C. Fort, F. Minardi, M. Inguscio, M. L. Chiofalo, and M. P. Tosi, Phys. Rev. Lett. 86, 4447 (2001).
• (9) F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999).
• (10) L. Pezzè, L. Pitaevskii, A. Smerzi, Stringari, G. Modugno, E. de Mirandes, F. Ferlaino, H. Ott, G. Roati, and M. Inguscio, Phys. Rev. Lett. 93, 120401 (2004).
• (11) C. Trallero-Giner, J. Drake, V. López-Richard, C. Trallero-Herrero, Joseph L. Birman, Physics Letters A 354, 115 (2006).
• (12) S. G. Mikhlin and K. L. Prössdorf, Approximate Methods for Solutions of Differential and Integral Equations (American Elsevier Publ. Co., NY, 1967).
• (13) I. G. Petrovskii, Lectures on the Theory of Integral Equations (Graylock Press, Rochester, 1957).
• (14) C. Trallero-Giner, J. Drake, V. López-Richard, C. Trallero-Herrero, Joseph L. Birman, Physica D, 237, 2342 (2008).
• (15) From general arguments it is possible to show that the expansion (4) converges in the mean, i.e. (see S. G. Mikhlin, Variational Methods in Mathematical Physics (Pergamon Press, 1964)). This allows to manipulate the integral in (3) and the series (4) such as to derive the relation (5).
• (16) H. Pu and N. P. Bigelow, Phys. Rev. Lett. 80, 1130 (1998).
• (17) H. Ott, E. de Mirandes, F. Ferlaino, G. Roati, G. Modugno, and M. Inguscio, Phys Rev. Lett. 92, 160601 (2004).
• (18) R. D. Lord, J. London Math. Soc. 24, 101 (1949).
• (19) I. S. Gradshteyn and I. M. Ryzhik, Tables of Integrals, Series and Products (Academic, NY, 1980).
• (20) Handbook of Mathematical Functions, edited by M. Abramowitz and I. Stegun (Dover, NY, 1972)
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters