Bose–Einstein condensates: Analytical methods for the Gross–Pitaevskii equation

Bose–Einstein condensates: Analytical methods for the Gross–Pitaevskii equation

Physics Letters A 354 (2006) 115–118 www.elsevier.com/locate/pla Bose–Einstein condensates: Analytical methods for the Gross–Pitaevskii equation Carl...

271KB Sizes 0 Downloads 15 Views

Physics Letters A 354 (2006) 115–118 www.elsevier.com/locate/pla

Bose–Einstein condensates: Analytical methods for the Gross–Pitaevskii equation Carlos Trallero-Giner a , J. Drake a , V. Lopez-Richard b,∗ , C. Trallero-Herrero c , Joseph L. Birman d a Faculty of Physics, Havana University, 10400 Havana, Cuba b F.F.C.L.RP, Departamento de Física e Matemática, Universidade de São Paulo, 14040-901 Ribeirão Preto, SP, Brazil c Physics and Astronomy Department, SUNY at Stony Brook, NY 11794-3800, USA d Department of Physics, The City College of CUNY, New York, NY 10031, USA

Received 30 June 2005; received in revised form 4 January 2006; accepted 11 January 2006 Available online 23 January 2006 Communicated by A.R. Bishop

Abstract We present simple analytical methods for solving the Gross–Pitaevskii equation (GPE) for the Bose–Einstein condensation (BEC) in the dilute atomic alkali gases. Using a soliton variational Ansatz we consider the effects of repulsive and attractive effective nonlinear interactions on the BEC ground state. We perform a comparative analysis of the solutions obtained by the variational Ansatz, the perturbation theory, the Thomas– Fermi approximation, and the Green function method with the numerical solution of the GPE finding universal ranges where these solutions can be used to predict properties of the condensates. Also, a generalization of the soliton variational approach for two-species of alkali atoms of the GPE is performed as a function of the effective interaction λi (i = 1, 2) and the inter-species λ12 and λ21 constants. © 2006 Elsevier B.V. All rights reserved.

Bose–Einstein condensates (BEC) have been reported in Li, Rb, and Na atoms [1]. In most cases the weakly interacting boson of the Gross–Pitaevskii (GPE) theory has been invoked to describe the properties of the condensate [2,3]. Nowadays, one of the most spectacular research is devoted to discuss the dynamics of BEC and the appearance of solitons based upon the nonlinear Schrödinger equation. The first bright soliton was created for negative scattering length BEC [4,5] and, for repulsive interaction, dark solitons have also been observed [6]. Recently, Eiermann et al. [7] reported experimental observation of different type of bright matter wave solitons for 87 Rb atoms in a periodic potential, allowing to study the dynamics of the condensate in a quasi-one-dimensional waveguide. Most of the theoretical work has been devoted to implement numerical solutions of the GP equation for the order parameter (see [3] and references therein). To design and to control the parameters of the condensate it is very useful to get analytical expressions for the chemical potential and for the order parameter as well. * Corresponding author.

E-mail address: [email protected] (V. Lopez-Richard). 0375-9601/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.physleta.2006.01.032

Our goal is to describe the dynamics of the dilute ultra-cold atom cloud in the BE condensed phase, implementing new analytical approaches beyond the Thomas–Fermi (TF) limit and perturbation theory [8]. In order to get analytical expressions, several authors have used a Gaussian as a trial wave function since in the linear limit one can obtain the linear Schrödinger equation in a parabolic potential [9]. For the negative scattering length case, a Gaussian function should probably be a good Ansatz to describe the condensate. Nevertheless, for repulsive atom–atom interaction, the shape of the wave function is similar to the TF solution (at least in the strong nonlinear limit). As it is well known, the validity of the variational result is in many cases qualitative and even if the shape of the order parameter is near to the trial wave function, this method is not always a good procedure for solving a nonlinear equation. In this Letter we present an efficient algorithm for the calculation of the ground state of the nonlinear Schrödinger equation (NLSE) in a harmonic trap using the Green function formalism and also a “soliton” variational solution. For sake of simplicity and in order to validate methodologically the obtained solutions, we just consider the isomorphic one-dimensional nonlinear GPE [2,3],

116

C. Trallero-Giner et al. / Physics Letters A 354 (2006) 115–118

describing the order parameter Φ, which can be written as h¯ 2 d 2 Φ 1 + mω2 x 2 Φ + λ|Φ|2 Φ = μΦ (1) 2m dx 2 2  with the normalization condition 1 = dx |Φ|2 . Here, μ is the chemical potential, ω is the trap frequency [10], and λ is an effective one-dimensional parameter that takes into account the interaction between the particles with mass m. In few cases, Eq. (1) allows an analytical solution. For example: (i) If the interaction λ is weak in comparison with the typical harmonic oscillator energy, the solution can be reached by perturbation in terms of the oscillator wave functions. (ii) If the kinetic energy can be neglected, in comparison with the self-interaction term, the Thomas–Fermi approximation applies [11]. (iii) For large negative self-interaction term, in comparison with the energy trap, the GPE yields the NLSE without an external potential, which presents a soliton solution (see Eq. (2)). A general method valid for λ > 0 and λ < 0 could provide a simple picture of the inside of the BEC. To get an explicit analytical solution of Eq. (1) for arbitrary values of λ we propose a variational Ansatz function of the type  1/2 K sech(Kx), ΦS (x, K) = (2) 2 −

where K is a variational parameter. The variational principle provides for the chemical potential μ[K] the expression 2 h¯ 2 K 2 π 2 mω2 Kλ + + , 3 2m 3 2K 2 3 where K must fulfill the transcendental equation μS [K] =

(3)

b4 − b3 − δ = 0

(4)

√ with b δ and l = h¯ /mω. The ¯ ¯ magnitude δ measures the ratio between the trap frequency ω and the coupling parameter λ in the condensate. Eq. (4) is a fourth-order algebraic polynomial problem with only one physically meaningful root. Using the Green function formalism, the GPE can be solved. Formally, Eq. (1) can be replaced by an integral equation [12] = K h2 /(mλ),

∞ Φ(x) =

= (l hω/λ)4 π 2 /4

 2  μG(x, x  )Φ(x  ) − λΦ(x  ) Φ(x  ) dx  ,



(5)

−∞

where the symmetric kernel G(x, x  ) is the Green function for the harmonic oscillator operator. It can be shown that the order parameter Φ(x) can be expressed in terms of the oscillator wave function {ϕn (x)} as Φ=



ϕn (x)Cn (μ),

(6)

n=0

where the vector coefficient solution C(μ) must fulfill the nonlinear system of equations

λ ¯ (μ) + (7) C ∗ T ∗ C C = 0, h¯ ωl

Fig. 1. Chemical potential in units of h¯ ω as a function of dimensionless self-interaction parameter λ/(h¯ ωl). Solid line: Green function formalism (7). Dashed line: Soliton variational approach (2). Dot: Thomas–Fermi approximation. Dash–dot: Perturbation theory. Empty circles: Numerical solution.

nm = (n + 1/2 − μ/h¯ ω)δnm (n = 0, 1, . . .) and T is a fourth order tensor [13]. The order parameter Φ and the chemical potential μ as a function of the dimensionless parameter λ/(h¯ ωl) can be obtained by solving the above exact nonlinear equation system (7). We have calculated the ground state solution Φ0 in terms of a truncated basis set, {ϕn (x)} (n = 1, . . . , I ), defining a finite-dimensional nonlinear Hill determinant eigenvalue equation [14] (I )  M μ[λ], C = 0, (I )

(8)

where Mnm (n, m = 1, 2, . . . , I ) are the corresponding matrix elements defined accordingly to Eq. (7). Much care must be taken in solving Eq. (8) because of the nonlinearity term ¯ ∗ T ∗ C. A judicious use of the general mathematical propC erties of the tensor T allows to get an efficient algorithm of resolution of Eq. (7) [13]. To solve Eq. (8), we implemented the von Neumann iterative procedure in a finite basis of di(0) mension I . The starting coefficients Cnm were chosen equal to δnm . For a particular λ value, the desired accuracy of μ defines (I ) the order of the matrix Mij . For positive values of λ < 10, a 25 × 25 matrix was enough to obtain μ/h¯ ω with a relative accuracy  = 10−6 . However, for negative values of λ, we obtain a lower convergence rate and a 50 × 50 matrix must be taken in order to obtain similar relative accuracy for λ > −10. Fig. 1 shows the dependence of the chemical potential on λ. All the physical parameters have been embedded in the dimensionless μ/h¯ ω and λ/(h¯ ωl). The result, following the Thomas–Fermi approach (indicated by dots), is compared with the solutions obtained by: (i) the soliton variational method (shown by dashed lines), (ii) perturbation theory (represented by short dasheddots), (iii) the exact solution using the Green function formalism (solid line). It can be seen that the TF approaches to, in the strong repulsive limit (λ > 0), the variational calculation and

C. Trallero-Giner et al. / Physics Letters A 354 (2006) 115–118

117

Fig. 2. μ/h¯ ω for large values of |λ/(h¯ ωl)| parameter. (a) Positive selfinteraction. (b) Negative self-interaction. Solid line: Green function formalism. Dashed line: Soliton variational approach. Dot: Thomas–Fermi approximation. Empty circles: Numerical solution.

the exact Green function solution. Fig. 2 shows the dependence of μ for large values of the self-interaction parameter in the condensate. We have compared our solutions with that obtained by direct numerical integration of Eq. (1). For the numerical solution we used a finite difference approximation transforming the GPE in an eigenvalue problem for a tri-diagonal matrix (for details see Ref. [15]). In Figs. 1 and 2 the numerical solutions are indicated by empty dots. From the graphics we extracted the following information: (i) The perturbation theory is remarkably good in the interval −1 < λ/(h¯ ωl) < 1.5. (ii) For λ = 0 the exact solution is μ/h¯ ω = 0.5, whereas our variational solution yields μ/h¯ ω = π/6  0.52, giving a difference of about 4%. (iii) For large repulsive interaction case [Fig. 2(a)], TF, Green function formalism, and soliton solutions agree very well. (iv) In the universal range of values −10 < λ < 10 and in comparison with the numerical calculations, we observe that the Green function formalism presents the accurate solution. (v) TF approaches asymptotically the numerical solution and Green function formalism for large values of the self-interaction term, while the variational solution has a discrepancy less that 4%. Also, we can see that the numerical solution and the Green function method match very well in the considered universal range of values of λ/(h¯ ωl). The variational soliton model presents an inexpensive estimation of the properties of the condensate in the range of self-interaction parameter −3 < λ/(h¯ ωl) < 10, where we have provided an analytical solution given by Eqs. (7) and (8) and its numerical implementation is readily simple. In order to analyze the quality of the order parameter Φ as a function of nonlinearity, we have evaluated the function  η(λ) = |Φnum −Φ|2 dx, where Φnum is the numerical solution of Eq. (1). The function η(λ) gives an estimation of the total error we have in the whole interval −∞ < x < ∞. For negative values of λ, the soliton variational approach presents, with respect to the numerical solution, a maximum integral difference of the order η(λ) = 0.3, while for λ > 0 this difference grows to 0.5. In the case of the Thomas–Fermi method, η(λ) ≈ 0.3 for

Fig. 3. Density profile l|Φ0 (x/ l)|2 as a function of the dimensionless selfinteraction parameter λ/(h¯ ωl). Upper panel: Variational solution. Lower panel: Thomas–Fermi approach.

λ/(h¯ ωl) = 5 decreasing as λ increases. If we consider the order parameter obtained by the Green function method and compare it with Φnum we obtain an η(λ) < 0.01 for the attractive interaction, while a maximum value of 0.035 is reached in the repulsive case. Fig. 3 displays a 3D picture of the single particle density profile l|Φ0 (x/ l)|2 according to the soliton variational approach (λ < 0 and λ > 0) and TF limit (λ > 0). We should stress that the plots represent a universal picture on the condensate order √ parameter lΦ0 (x) in terms of the dimensionless parameter λ/(h¯ ωl). It can be seen that the wave function is more delocalized and the maximum of l|Φ0 (x/ l)|2 decrease as λ increases for positive values, while the opposite is observed for the attractive interaction case. The above results can be straightforward generalized to the case of two-species in the condensate. Here, we have to solve two coupled GPE where μi , ωi , mi , Φi , and λi are the parameter for the ith species. The terms describing the interaction between both species are given by λ12 and λ21 and can be positive or negative magnitudes. We just considered the variational procedure. The trial normalized wave functions of the soliton types are  ΦiS (x, Ki ) =

Ki 2

1/2 sech(Ki x)

(9)

118

C. Trallero-Giner et al. / Physics Letters A 354 (2006) 115–118

Fig. 4. Two-species. μ1 /h¯ ω1 as a function of λ1 /(h¯ ω1 l1 ) for several values of the dimensionless two-species interaction η = λ21 /(l2 hω ¯ 2 ). In the calculation λ2 /(h¯ ω2 l2 ) = 7.5, l2 /1l = 3/2.

containing two unknown variational parameters Ki (i = 1, 2). The required values of Ki correspond to the minimum of μi [K1 , K2 ], leading to the equations    3 bi 3 λ¯ ij λ¯ i li 4 3 F − δi = 0, i = j b i + bi 1 + 4 λ¯ i λ¯ j lj bj with bi = Ki h¯ 2 /(mi λi ), δi = ( π2 λ¯1 )4 , and i ∞ F (Z) =

λ¯ i = λi /(li h¯ ωi ), λ¯ ij = λij /(li h¯ ωi ),



2  sech(x) sech(Zx) 1 − Zx tanh(Zx) dx.

−∞

For the chemical potentials we obtain   3  

 bi 1 λ¯ i li 3 λ¯ ij 1 3δi μi = − λ¯ 2i bi2 1 − H − 2 , 6 2 λ¯ i bi h¯ ωi λ¯ j lj bj bi where ∞ H (Z) =

 2 Zx sech(x) sech(Zx) tanh(Zx) dx.

−∞

Fig. 4 shows the dependence of the chemical potential μ1 , in units of the energy trap h¯ ω1 , on the dimensionless selfinteraction parameter λ1 /(h¯ ω1 l1 ) for, λ2 /(h¯ ω2 l2 ) = 7.5 and several values of the ratio λ12 /λ21 . In the calculation we chose 5 λ21 λ12 = l1 hω ¯ 1 2 l2 h¯ ω2 and l2 / l1 = 3/2. From Fig. 4, one can see the peculiar effects of the interaction between the two-species on μ1 . As the interaction λ12 becomes negative, μ1 decreases and the condensate becomes more confined. In the opposite case, when λ12 > 0 if

λ1 is positive, μ1 increases; while for λ1 < 0, the chemical potential μ1 is practically independent of the interaction between both species. In summary, the reported analytical techniques allow us to explore regions of positive, negative, and zero nonlinear interaction in condensates. We have been able to estimate the range of applicability of the perturbation theory, TF, soliton variational calculation and Green function formalism through a unique universal parameter λ/(h¯ ωl) (see Figs. 1 and 2). These can be useful tools to design different traps for one and twospecies of condensates. Since we obtained compact solutions in a soliton like manner, we are also able to the study the dynamics of the BEC. The ground state wave functions of two-species of BEC have been obtained in the past within the TF approximation [16] and by numerical methods [15]. We, for the first time, provided a simple analytical solution for the order parameters and the chemical potentials of double condensates valid for any values of effective interaction constants λi (i = 1, 2) and effective interactions λ12 and λ21 between both species. The generalization of the methods here reported to realistic 3D problems is straightforward and will allow us to predict quantitative phenomena, to investigate the formation of vortices or to describe the collisional relaxation times or self compression, among other novel properties of the binary mixtures. Acknowledgements This work was supported in part from the Science Division of the The City College of CUNY and from the CUNYCaribbean Exchange Program. V.L.R. acknowledges FAPESP (Brazil) for the financial support. References [1] C.C. Bradley, C.A. Sackett, J.J. Tollett, R.G. Hulet, Phys. Rev. Lett. 75 (1995) 1687; J. Han, et al., Phys. Rev. A 57 (1998) R4114; K.B. Davis, et al., Phys. Rev. Lett. 75 (1995) 3969. [2] V.L. Ginzburg, L.P. Pitaevskii, Zh. Eksp. Teor. Fiz. 34 (1958) 1240; E.P. Gross, J. Math. Phys. 4 (1963) 195. [3] F. Dalfovo, S. Giorgini, L.P. Pitaevskii, S. Stringari, Rev. Mod. Phys. 71 (1999) 463. [4] L. Khaykovich, et al., Science 296 (2002) 1290. [5] K.E. Strecker, et al., Nature (London) 417 (2002) 150. [6] S. Burger, et al., Phys. Rev. Lett. 83 (1999) 5198. [7] B. Eiermann, et al., Phys. Rev. Lett. 92 (2004) 230401. [8] J.-R. Yan, J. Lung, S.-M. Ao, D.-B. Cao, Chin. Phys. Lett. 19 (2002) 1245. [9] V.M. Pérez-García, et al., Phys. Rev. A 56 (1997) 1424. [10] M.H. Anderson, et al., Science 269 (1995) 198. [11] M. Edwards, K. Burnett, Phys. Rev. A 51 (1995) 1382; G. Baym, C. Pethick, Phys. Rev. Lett. 76 (1996) 6. [12] P.M. Morse, H. Fesbach, Methods of Theoretical Physics, McGraw–Hill, New York, 1953. [13] A detailed analysis of the properties and symmetry of the matrix elements Tp,l,n,m are given in C. Trallero-Giner, et al., unpublished. [14] C.M. Bender, S.A. Orszag, Advanced Mathematical Methods for Scientists and Engineers, McGraw–Hill, New York, 1978. [15] H. Pu, N.P. Biegelou, Phys. Rev. Lett. 80 (1998) 1130. [16] T.-L. Ho, V.B. Sheniy, Phys. Rev. Lett. 77 (1996) 3276.