Spectral method for the time-dependent Gross-Pitaevskii equation with a harmonic trap

Spectral method for the time-dependent Gross-Pitaevskii equation with a harmonic trap

Claude M. Dion and Eric Cance`s

CERMICS, E´ cole Nationale des Ponts et Chausse´es, 6 & 8 avenue Blaise Pascal, Cite´ Descartes, Champs-sur-Marne, 77455 Marne-la-Valle´e, France

共Received 5 July 2002; published 15 April 2003兲

We study the numerical resolution of the time-dependent Gross-Pitaevskii equation, a nonlinear Schro¨dinger equation used to simulate the dynamics of Bose-Einstein condensates. Considering condensates trapped in harmonic potentials, we present an efficient algorithm by making use of a spectral-Galerkin method, using a basis set of harmonic-oscillator functions, and the Gauss-Hermite quadrature. We apply this algorithm to the simulation of condensate breathing and scissor modes.

DOI: 10.1103/PhysRevE.67.046706 PACS number共s兲: 02.70.Hm, 03.75.Kk, 31.15.⫺p I. INTRODUCTION

The experimental realization of Bose-Einstein condensa- tion 关1–3兴 has prompted much work on the study of the dynamics of these condensates. From the theoretical side, many interesting results have been obtained using the Gross- Pitaevskii equation共GPE兲 关4–6兴,

iប⳵⌿

⳵^{t ⫽}冋^{⫺2mប2 ⵜ2⫹Vext⫹4␲បm2aN兩⌿兩2}册^{⌿, 共1兲}

with the normalization condition 储⌿(t)储L²⫽1᭙ t, to de- scribe the order parameter ⌿ 共also called the condensate wave function兲 of N condensed bosons of mass m, interacting via a contact potential described by the scattering length a, and eventually confined by an external potential V_ext. Even though the Gross-Pitaevskii equation is based on the ap- proximation that all bosons are in the condensed phase (T

⫽0 K), direct comparison between theoretical and experi- mental results have shown that, in many cases, solutions of the GPE contain the essential physics of the underlying phe- nomena 关7–10兴. This nonlinear Schro¨dinger equation 共NLSE兲 has been used, in its time-dependent form, to inves- tigate many aspects of the dynamics of Bose-condensed gas, such as the formation of vortices 关11兴, the interference be- tween condensates 关12兴, of the possibility of creating atom lasers关13,14兴, to mention only a few.

Most of these and other numerical studies of the time- dependent GPE are based on grid methods, i.e., discretize the spatial coordinates on a grid of points, the resulting differen- tial equation being usually solved by Crank-Nicholson or split-operator Fourier methods共see, e.g., Refs. 关15–18兴兲. We must point out that, while much care must be taken in solv- ing Eq. 共1兲 because of the nonlinearity, we find, to our dis- may, that many authors give results calculated with the time- dependent GPE without even specifying what method they have used for their numerical simulation.

In this paper, we wish to focus our attention on the case where the Bose-Einstein condensate is in a 共possibly aniso- tropic兲 harmonic trap, i.e.,

V_ext共X,Y,Z兲⫽1 2m共␻x

2X²⫹␻y 2Y²⫹␻z

2Z²兲, 共2兲

which is the case for most experimental setups关19,20兴. The method we propose is based on the spectral decomposition of

⌿ on a basis of harmonic-oscillator wave functions. In such a representation, the kinetic ⫹ trapping potential part of the Hamiltonian is diagonal. The nonlinear part is computed by forward and backward transformations from the spectral to a grid representation. By judicious use of the Gauss-Hermite quadrature, this can lead to an algorithm that is more effi- cient than those based on grid methods. Although this is akin to discrete variable representation 共DVR兲 methods based on Hermite polynomials, which have been successfully used for the time-independent and time-dependent GPE 关21,22兴, our method is distinct, since our Hamiltonian is expressed in the spectral representation for both the kinetic and potential op- erators.

We expose in Sec. II our spectral method and the resulting algorithm. We then present different time-evolution schemes that can be used in combination with the spectral method. We finally give in Sec. IV some results that can be obtained from the numerical simulation of the time-dependent GPE, namely, the study of condensate breathing and scissor modes.

II. SPACE DISCRETIZATION

To simplify the calculation, we will first rescale Eq.共1兲 in the three spatial dimensions (X,Y ,Z) and in time,

X⫽冉^mប␻x冊^{1/2x, 共3a兲}

Y⫽冉^mប␻y冊^{1/2y , 共3b兲}

Z⫽冉^mប␻z冊^{1/2z, 共3c兲}

t⫽ 1

␻z␶^. 共3d兲

We also introduce a new wave function ␺ ^{defined as}

⌿共t,X,Y,Z兲⫽A␺共␶^{,x, y ,z}兲, and, considering the normalization condition

冕_R³兩⌿共t,X,Y,Z兲兩²dXdY dZ⫽1᭙t, we choose

A⫽冉^mប冊^{3/4共␻x␻y␻z兲1/4,}

such that

冕R³兩␺共␶^{,x,y ,z}兲兩²dxdy dz⫽1.

The Gross-Pitaevskii equation therefore becomes

i⳵␺

⳵␶^⫽冋^␻␻xz冉^{⫺12ⵜx2⫹x22}冊^⫹␻␻yz冉^{⫺12ⵜy2⫹y22}冊

⫹冉^{⫺12ⵜz2⫹z22}冊^{⫹␭兩␺兩2}册^{␺, 共4兲}

with

␭⫽4␲^aN冉^{mប ␻␻x␻zy}冊^{1/2. 共5兲}

Coordinate z should be chosen such that␻zis the greatest of the three frequencies关this is related to the arbitrary choice of the scaling factor in Eq.共3d兲兴.

As all the physical parameters have been absorbed in the nonlinear parameter␭, calculations with the same ␭ can cor- respond to results for different species, but in diverse experi- mental conditions. We can define acceptable lower and upper bounds for ␭ by considering the effective range of the dif- ferent physical parameters. Considering only cases where the interparticle interaction is repulsive, i.e., a⬎0 and therefore

␭⬎0, at the lower end we can consider a small ⁴He* con- densate (m⫽4.0 amu, a⫽302 a.u. 关23兴兲 of N⫽10³ atoms in a highly anisotropic ␻x␻y/␻z⫽2␲⫻10^⫺1Hz trap, giving

␭⬇1.3, while for a bigger N⬃10⁶ condensate of heavy at- oms such as ⁸⁷Rb (m⫽86.9 amu, a⫽106 a.u. 关24兴兲, ␭ can reach 10⁵ for isotropic traps. In the following, we will re- strict our study to␭ in the range 1 –10³, considering that the Thomas-Fermi approximation can be used for greater values of ␭ 关21兴.

A. The spectral-Galerkin method in 1D

For pedagogical purposes, we first explain our numerical method on the simple case of the one-dimensional 共1D兲 NLSE

i⳵␺

⳵^{t 共t,x兲⫽H0}␺共t,x兲⫹␭兩␺共t,x兲兩²␺共t,x兲, 共6兲 with

H₀⫽⫺1 2

⳵²

⳵^x2⫹ 1 2x².

Extensions to the three-dimensional case will be detailed in the following section.

Denoting by ␺(t) the function x哫␺(t,x), it can be proven关25兴 that if

␺0苸Xª再^{␹苸L2共R兲, 冕R冏⳵␹⳵x冏2⬍⫹⬁,}

冕_R^{x2兩␹共x兲兩2dx⬍⫹⬁}冎^,

Eq. 共6兲 with initial condition ␺0 has a unique solution in C⁰(关0,⫹⬁关,X)艚C¹„关0,⫹⬁关,L²(R)… and that both the L² norm

储␺共t兲储L²⫽冋^{冕R兩␺共t,x兲兩2dx}册^1/2,

and the energy

E⫽„H0␺共t兲,␺共t兲…⫹␭

2冕R兩␺共t,x兲兩⁴dx

are conserved by the dynamics. A variational formulation of Eq. 共6兲, supplemented by the initial condition ␺^(t⫽0)

⫽␺0, where␺0苸X reads

Search ␺苸C⁰共关0,T兴,X兲艚C¹„关0,T兴,L²共R兲… such that ᭙␹苸X, i d

dt „␺共t兲,␹…⫽„H0␺共t兲,␹…⫹␭„兩␺共t兲兩²␺共t兲,␹…, 共7兲

␺共0兲⫽␺0.

Numerical solutions can then be obtained by approximating problem 共7兲 with a Galerkin method: a finite-dimensional subspaceXNof the infinite-dimensional vector spaceX being given, we consider

Search ␺N苸C¹共关0,T兴,XN兲 such that ᭙␹N苸XN,

i d

dt „␺N共t兲,␹N…⫽„H0␺N共t兲,␹N…⫹␭„兩␺N共t兲兩²␺N共t兲,␹N…, 共8兲

␺N共0兲⫽␺0.

Denoting by (␾0, . . . ,␾N) an orthonormal basis of X_N for the L² scalar product and by C(t)⫽关cn(t)兴0⭐n⭐N the vector of C^N⫹1 collecting the coefficients of ␺N(t) in the basis (␾0, . . . ,␾N), i.e.,

␺N共t,x兲⫽_n兺_⫽0^{N cn共t兲␾n共x兲,}

problem共8兲 can be reformulated as

Search C苸C¹共关0,T兴,C^N⫹1兲 such that

共2003兲

idC

dt 共t兲⫽hC共t兲⫹␭F„C共t兲…, C共0兲⫽C0, 共9兲

where C₀are the coefficients of␺0and h the matrix of H₀in the basis (␾0, . . . ,␾N)

关C0兴n⫽共␺0,␾n兲L², h_nm⫽共H0␾m,␾n兲, and where the function F is defined by

F共C兲n⫽_k,l,m兺^N_⫽0^{Iklmnck*clcm, 共10兲}

with

I_klmn⫽冕R␾k*␾l␾m␾n*.

The efficiency of a direct implementation 关26,27兴 of the Galerkin method described above is very poor: the calcula- tion of the integrals I_klmn 共which can be precomputed if the basis is small enough that the integrals can be stored in memory兲 scales as O(N⁴N_p), where N_p is the number of grid points of the quadrature method, and the computation cost for one evaluation of the function F scales as N⁴ 关for each of the N coefficients, O(N³) operations are needed兴.

Our aim is to show that the Galerkin method becomes very efficient if (␾0, . . . ,␾N) are the N⫹1 lowest eigen- modes of the harmonic oscillator H₀. In this case, indeed, the vector F(C) can be computed exactly 共up to round-off er- rors兲 in O(N²) operations. Let us recall that the eigenmodes (␾n)_n苸Nof H₀ read

␾n共x兲⫽Hn共x兲e^⫺x2/2,

where Hn(x) is the nth Hermite polynomial 关28兴, and that they satisfy

H₀␾n⫽En␾n with E_n⫽n⫹1 2.

In such a basis, the matrix h is therefore diagonal: h

⫽diag(E0, . . . ,E_n). In addition, for any C苸C^N⫹1, one has

F共C兲n⫽冕R兩␺共x兲兩²␺共x兲␾n共x兲dx, 共11兲 where ␺^(x)⫽兺nN⫽0

c_n␾n(x). The key point is now that for any n⭐N the integrand in Eq. 共11兲 is of the form Q(x)e^⫺2x2, where Q(x) is a polynomial of degree lower or equal to 4N; each of the N⫹1 integrals can therefore be computed exactly with a Gauss-Hermite quadrature formula involving 2N Gauss points 关29兴. More precisely, we have, for any polynomial Q of degree lower or equal to 4N,

冕⫺⬁

⫹⬁

Q共x兲e^⫺x2dx⫽^2N_k兺_⫽1^{⫹1 wkQ共xk兲,}

where兵^xk其 are the roots of the Hermite polynomial H2N⫹1

and where兵^wk其 are convenient weights关30兴. By a change of variable in integral共11兲, it follows that

F共C兲n⫽^2N_k兺_⫽1^⫹1 冉^w冑ke2xk2冊^{兩␺共xk/冑2兲兩2␺共xk/冑2兲␾n共xk/冑2兲.}

Spectral-Galerkin methods are usually not very efficient 关31兴; but they can be in the specific case of the NLSE, we are interested in because of the special form of the nonlinearity.

Let us now denote by P苸M(N⫹1,2N⫹1) the matrix collecting the values of the basis functions (␾n)_0⭐n⭐Nat the Gauss points (x_k)_{1⭐k⭐2N⫹1}:

P_nk⫽␾n共xk/冑²兲,

and by w˜_k⫽wke^xk2/冑2. An efficient algorithm for the com- putation of F(C) for a given C苸C^N⫹1 reads the following.

共1兲 Compute the vector ⌿苸C^2N⫹1 defined by

⌿⫽P^T•C.

共2兲 Compute the vector ⌶苸C^2N⫹1 coefficient by coeffi- cient along formula

⌶k⫽w˜k兩⌿k兩²⌿k. 共3兲 Compute

F共C兲⫽P•⌶.

The vectors C and ⌿ are the representation of the wave function␺in the spectral basis兵␾n其0⭐n⭐N and in real space 共at the 2N⫹1 Gauss points兵^xk/冑²其), respectively. Steps共1兲 and共3兲 of the above algorithm scale quadratically in N 共these are matrix-vector products兲, and step 共2兲 scales linearly in N.

We therefore end up with an algorithmic complexity in O(N²).

In practice, the function C哫F(C) is called one or several times at each time step; of course, the matrix P as well as the weights w˜

k can be precomputed once and for all and stored in memory.

B. The spectral-Galerkin method in 3D

Let us now turn to the 3D setting and consider the res- caled equation

i⳵␺

⳵^{t 共t,x,y,z兲⫽}冋^{␻␻xzH0共x兲⫹␻␻yzH0共y兲⫹H0共z兲}册^{␺共t,x,y,z兲}

⫹␭兩␺共t,x,y,z兲兩²␺共t,x,y,z兲, 共12兲 with

H₀共x兲⫽⫺1 2

⳵²

⳵^x2⫹ 1

2x², H₀共y兲⫽⫺1 2

⳵²

⳵^y2⫹ 1 2y²,

H₀共z兲⫽⫺1 2

⳵²

⳵^z2⫹ 1 2z².

For␭⭓0, a global-in-time existence and uniqueness result is available for Eq.共12兲 with initial condition␺^(t⫽0)⫽␺0and

␺0苸X⫽兵␹^苸L2共R³兲, “␹苸关L²共R³兲兴³, 共x²⫹y²⫹z²兲^1/2␹苸L²共R³兲其^.

On the other hand, it is well known that finite-time blow up may be observed for ␭⬍0 and for some initial conditions 关25兴. As stated above, we focus here on the case where ␭

⭓0.

Following the same lines as in the Sec. II A, the approxi- mated wave function␺N(t) is expended on the spectral ten- sor basis set

„␾n_x共x兲␾n_y共y兲␾n_z共z兲…0⭐n_x⭐N_x,0⭐n_y⭐N_y,0⭐n_z⭐N_z.

One therefore has

␺N共t,x,y,z兲

⫽n兺_x⫽0 N_x

n兺_y⫽0 N_y

n兺_z⫽0 N_z

c_n

xn_yn_z共t兲␾n_x共x兲␾n_y共y兲␾n_z共z兲.

共13兲 The equation satisfied by the three index tensor C

⫽关cn_xn_yn_z兴 in the Galerkin approximation formally has the same expression as in 1D,

idC

dt 共t兲⫽hC共t兲⫹␭F„C共t兲…, the linear operator h now being defined by

关hC兴n_xn_yn_z⫽En_xn_yn_zc_n

xn_yn_z, with

E_n

xn_yn_z⫽␻x

␻z冉^nx⫹12冊^⫹␻␻yz冉^ny⫹12冊^⫹冉^nz⫹12冊^,

and the nonlinear function F(C) by

F共C兲]n_xn_yn_z⫽冕R³兩␺共x,y,z兲兩²␺共x,y,z兲␾n_x共x兲␾n_y共y兲

⫻␾n_z共z兲dxdydz,

where␺(x, y ,z) is given by Eq.共13兲.

Let us denote by 兵^xk其¹⭐k⭐2Nx⫹1, 兵^yk其¹⭐k⭐2Ny⫹1, 兵^zk其¹⭐k⭐2Nz⫹1 the roots of the Hermite polynomials H2N_x⫹1, H2N_y⫹1, H2N_z⫹1 and 兵^wk

x其¹⭐k⭐2N_x⫹1, 兵^wk

y其1⭐k⭐2N_y⫹1, 兵^wk

z其1⭐k⭐2N_z⫹1 the associated summation weights. Let us also introduce the matrices P_x苸M(Nx

⫹1,2Nx⫹1), Py苸M(Ny⫹1,2Ny⫹1), Pz苸M(Nz⫹1,2Nz

⫹1) defined by

关Px兴n_xk_x⫽␾n_x共xk_x/冑²兲, 关Py兴n_yk_y⫽␾n_y共yk_y/冑²兲, 关Pz兴n_zk_z⫽␾n_z共zk_z/冑²兲,

and the weights

w˜

k_x x⫽w_k

x xe^xkx

2

冑^{2 , w˜kyz⫽}

w_k

y ye^yky

2

冑^{2 , w˜kzz⫽}

w_k

z z e^zkz

2

冑^{2 .}

The following algorithm for the computation of F(C) scales in O(NN_xN_yN_z) where N⫽max(Nx,N_y,N_z).

共1兲 Set ⌿^SSS⫽C.

共2兲 Compute ⌿n_xn_yk_z SSR ⫽n兺_z⫽0

N_z

关Pz兴n_zk_z⌿n_xn_yn_z SSS ,

O共NxN_yN_z²兲 operations.

共3兲 Compute ⌿n_xk_yk_z SRR ⫽_n兺

y⫽0 N_y

关Py兴n_yk_y⌿n_xn_yk_z SSR ,

O共NxN²_yN_z兲 operations.

共4兲 Compute ⌿k_xk_yk_z RRR ⫽_n兺

x⫽0 N_x

关Px兴n_xk_x⌿n_xk_yk_z SRR ,

O共Nx

2N_yN_z兲operations.

共5兲 Compute ⌶k_xk_yk_z RRR ⫽w˜k_x

xw˜

k_y y w˜

k_z z兩⌿k_xk_yk_z

RRR 兩²⌿k_xk_yk_z RRR ,

O共NxN_yN_z兲 operations.

共6兲 Compute ⌶k_xk_yn_z RRS ⫽_k兺

z⫽1 2N_z⫹1

关Pz兴n_zk_z⌶k_xk_yk_z RRR ,

O共NxN_yN_z²兲 operations.

共7兲 Compute ⌶k_xn_yn_z RSS ⫽ k兺_y⫽1

2N_y⫹1

关Py兴n_yk_y⌶k_xk_yn_z RRS ,

O共NxN²_yN_z兲 operations.

共8兲 Compute ⌶n_xn_yn_z SSS ⫽ _k兺

x⫽1 2N_x⫹1

关Px兴n_xk_x⌶k_xn_yn_z RSS ,

O共Nx

2N_yN_z兲 operations.

共9兲 Set F共C兲⫽⌶^SSS.

In the above formulation, the superscripts S and R stand for spectral and real space representations, respectively. In other words, steps共2兲–共4兲 constitute the successive transformation of the wave function from the spectral basis to a spatial rep- resentation on the series of points of the Gauss-Hermite 共2003兲

quadrature. The nonlinear term of the Hamiltonian is then calculated in this spatial representation关step 共5兲兴, while steps 共6兲–共8兲 correspond to the inverse transform back to the spec- tral basis. It is this procedure of forward and backward trans- formation that allows us to obtain a much better scaling than the implementation of Eq.共10兲.

The scaling of the above algorithm 关O(N⁴) if N_x⫽Ny

⫽Nz] has to be compared with the scaling of fast Fourier transform based algorithms which scale in O„Np

3log₂(N_p)…, where N_p is the number of grid points per direction. The main interest of the spectral method is that for a similar accuracy, the number of spectral basis functions per direction 共here denoted by N) can usually be chosen much smaller than the number N_p of grid points per direction. This is es- pecially true when the problem considered displays a sym- metry in one or more of the directions, in which case the basis set used in the Galerkin approximation Eq.共13兲 can be restricted to even harmonic-oscillator functions共in the corre- sponding direction兲. We will come back on this important feature of the spectral method in Sec. IV.

C. Exploiting spherical or cylindrical symmetry When ␻x⫽␻y⫽␻z the one-particle Hamiltonian pos- sesses spherical symmetry. If the initial condition ␺0⫽␺^(t

⫽0) has the same symmetry, then the wave function␺^{(t) is} spherical symmetric for any t⬎0: ␺(t,x, y ,z)⫽␺^(t,r), where r⫽(x²⫹y²⫹z²)^1/2 is the radial coordinate. Equation 共4兲 leads to the effective 1D dynamics

i⳵␺

⳵^{t ⫽}冋^{⫺2r12 ⳵⳵r冉r2⳵⳵r冊⫹r22⫹␭兩␺兩2}册^{␺. 共14兲}

Let us now define the function

␹共t,r兲⫽再^{冑⫺2冑␲2r␲␺r共t,r兲␺}共t,⫺r兲 if r⬍0.^{if r⬎0} It is easy to check that␹ actually satisfies the 1D NLSE

i⳵␹

⳵^{t ⫽H0}␹⫹␭ 兩␹兩² 2␲^r2␹^.

Besides, for any t⬎0 the function ␹^{(t): r}哫␹(t,r) is odd and belongs to L²(R) since

冕⫺⬁

⫹⬁兩␹共t,r兲兩²dr⫽冕^{0⫹⬁4␲r2兩␺共t,r兲兩2dr⫽1.}

It can thus be expanded on the odd modes of the harmonic oscillator:

␹共t,r兲⫽_n兺^⫹⬁_⫽0 ^{cn共t兲␾2n⫹1共r兲.}

A spectral-Galerkin approximation can now be used. The vector C(t)苸C^N⫹1 collecting the coefficients „c^k(t)…⁰⭐k⭐N

of the approximated wave function

␹N共t,r兲⫽n兺⫽0 N

c_n共t兲␾2n⫹1共r兲,

obeys once again a dynamics of the form

idC

dt 共t兲⫽hC共t兲⫹␭F„C共t兲….

Here

h⫽diag共E2n⫹1兲 with E2n⫹1⫽2n⫹3 2 and

关F共C兲兴n⫽冕R

兩␹共r兲兩²

2␲^r2 ␹共r兲␾2n⫹1共r兲dr, where ␹^(r)⫽兺nN⫽0

c_n␾2n⫹1(r). As for any 0⭐n⭐N,

␾2n⫹1(r)⫽rP2n(r)e^⫺r2/2 where P_2n is a polynomial of de- gree equal to 2n, it follows that the above integrals can be computed exactly with 4N Gauss points.

Let us now turn to the cylindrical symmetry when 共for instance兲␻x⫽␻y and when the initial data reads␺0(x, y ,z)

⫽␺0(r,z) with r⫽(x²⫹y²)^1/2. In this case, the cylindrical symmetry is preserved by the dynamics so that for any t

⬎0, ␺(t,x, y ,z)⫽␺^(t,r,z) and the time evolution of

␺(t,r,z) is then governed by the 2D equation

i⳵␺

⳵^{t ⫽}冋^{␻␻xz再⫺2r1 ⳵⳵r冉r⳵⳵r冊⫹r22冎⫹}冉^{⫺12 ⳵⳵z22⫹z22}冊

⫹␭兩␺兩²册^{␺, 共15兲}

set on the spatial domain R^⫹⫻R. Defining a new function

␹(t,r,z) on the space domainR² by

␹共t,r,z兲⫽再^{␺␺共t,r,z兲}共t,⫺r,z兲 if r⬍0,^{if r⬎0} it occurs that␹ ^satisfies

i⳵␹

⳵^{t ⫽}冋^␻␻xz冉^{⫺12 ⳵⳵r22⫹r22}冊^⫹冉^{⫺12 ⳵⳵z22⫹z22}冊^{⫺␻␻xz 2r1 ⳵⳵r}

⫹␭兩␹兩²册^{␹ 共16兲}

on the space domainR², and that, by construction, the func- tion r哫␹(t,r,z) is even. A spectral-Galerkin approximation is obtained by expanding the wave function on the spectral tensor basis set

„␾2n_r共r兲␾n_z共z兲…0⭐n_r⭐N_r,0⭐n_z⭐N_z.

The coefficients (c_n

rn_z)_0⭐n

r⭐N_r,0⭐n_z⭐N_z of the expansion are solution of an equation of the same form as above,

idC

dt 共t兲⫽hC共t兲⫹␭F„C共t兲….

The main difference is that in this case, the linear map h takes into account the operator⫺(1/2r)(⳵^/⳵^r):

关hC兴n_rn_z⫽冋^␻␻xz冉^nr⫹12冊^⫹冉^nz⫹12冊册^Cnrnz

⫺1 2

␻x

␻z _m兺

r⫽0 N_r

冉^{1r d␾dr2mr,␾2nr}冊L²

C_m

rn_z.

Let us remark that the scalar product

„(1/r)(d␾2m_r/dr),␾n_r…L² is well defined since the first de- rivative of␾2n_ris of the form r P_2n

r(r)e^⫺r2/2where P_2n

ris a polynomial of degree 2n_r; in addition, it can be computed exactly by numerical integration with 2n_rGauss points. It is worth pointing out that the ‘‘Hamiltonian’’ in Eq.共16兲 is not self-adjoint because of the term ⫺(1/2r)(⳵^/⳵r) and that the L² norm of ␹(t) is not a conserved quantity; on the other hand, the L² norm of ␹(t) for the measure r dr dz is con- served.

III. TIME DISCRETIZATION

When a spectral-Galerkin method is used to discretize the space variables, one ends up with a finite-dimensional dy- namical system of the form

idC

dt 共t兲⫽hC共t兲⫹␭F„C共t兲…, 共17兲 with initial condition C(t⫽0)⫽C0. We then use a basic fourth-order Runge-Kutta method关32兴 to solve Eq. 共17兲. Let us mention that, as the Hamiltonian character of the NLSE is preserved by the spectral-Galerkin discretization, it would be possible to resort to symplectic methods 关33兴; such algo- rithms, which are particularly advised for long time evolu- tion, are, however, not tested in the present work.

We will also use a grid method, based on the split- operator method, to serve as a benchmark for the spectral algorithm we have just detailed. We recall below the main features of this approach.

The wave function at time ␶⫹⌬␶ can be obtained from the wave function at␶ according to

␺共␶⫹⌬␶兲⫽Uˆ共␶^,␶⫹⌬␶兲␺共␶兲, 共18兲 with the propagator Uˆ (␶^,␶⫹⌬␶) being expressed, for suffi- ciently small intervals⌬␶^{, as}

Uˆ共␶^,␶⫹⌬␶兲⫽exp关⫺iH共␶兲⌬␶兴, 共19兲 where H(␶) is the Hamiltonian of Eq.共4兲. As the potential and nonlinear components of the Gross-Pitaevskii Hamil- tonian do not commute with the kinetic operator, we apply the split-operator method关34兴 to obtain

exp关⫺iH共␶兲⌬␶兴⫽exp冋^⫺iT⌬2␶册^{exp关⫺i共V⫹␭兩␺兩2兲⌬␶兴}

⫻exp冋^⫺iT⌬2␶册^{⫹O共⌬␶3兲, 共20兲}

with T the kinetic operator and V the trapping potential. The middle term is diagonal in position space, while the kinetic part is diagonal in momentum space. A fast Fourier trans- form is thus used before application of the kinetic operator, followed by the inverse transform. Note that if the interme- diate wave function at time ␶⫹⌬␶ is not needed, the two successive kinetic operators half steps can be combined.

From a previous study关35兴, it appears that the split-operator method is the fastest algorithm for solving a NLSE on a grid.

IV. RESULTS

The first test we perform is the propagation of the ground stationary state 共obtained from the time-independent GPE solved by a method based on the optimal damping al- gorithm 关36–38兴兲, while monitoring the value of the coeffi- cients c(␶) of the expansion 共13兲. For the spherically sym- metric case, we require that the relative error on the c₀ coefficient 共which has the largest absolute value兲 be inferior to 10^⫺8, i.e.,兩兩c0(␶⁾兩²⫺兩c0(␶⫽0)兩²兩/兩c0(␶⫽0)兩²⭐10^⫺8᭙␶ 苸关0,100兴. This criterion also results in an absolute error of all coefficients 兩cn(␶⁾兩²⫺兩cn(␶⫽0)兩²⭐10^⫺8. We have also checked that the phase of the coefficients is correct, by cal- culating 兩cn(␶⁾⫺cn(␶⫽0)e^⫺i␮␶兩²/兩cn(␶⁾兩², where ␮ ^{is the} chemical potential of the ground stationary state of the GPE 关6兴, and this value indeed is less than 10^⫺12.

In this 1D case, we need N⫽20 basis functions for ␭

⫽100, and the resulting time step for the Runge-Kutta propagator is ⌬␶⫽0.005. If ␭⫽1000, the basis set used should be slightly larger, N⫽26, with a smaller time step

⌬␶⫽0.0025 to insure that the above error criteria are met.

The resulting propagation time up to ␶⫽100 is 8.9 s for ␭

⫽100 共calculated on an Athlon 1.2 GHz processor running under Linux, using the NAG Fortran 95 compiler at the⫺O2 level of optimization兲 and 28.3 s for ␭⫽1000. If we double the size of the basis set, we get a CPU time of 32.9 s for␭

⫽100, showing the expected O(N²) scaling of the algorithm in 1D.

Comparing now with the grid method described in Sec.

III, we use N_p⫽64 grid points in the range ⫺8⭐r⭐8. The time step used is ⌬␶⫽0.000 25, resulting in a propagation time of 10.3 s, which is slightly longer than what we obtain using the Runge-Kutta method.

We now apply our algorithm to study the dynamics of trapped condensates. Referring again to the spherically sym- metric case, we start with the stationary ground state for an isotropic trap frequency ␻. We then let this initial state␺0

evolve in a trap of frequency ␻/2, as illustrated in Fig. 1, corresponding to an experiment where the frequency of the potential trapping the condensate would be instantaneously reduced by a factor of 2. The corresponding time-evolving wave function 兩␺^(t,r)兩² is shown in Fig. 2, for ␭⫽10. We must note that the values of ␭ we give correspond to the 共2003兲

condensate in the initial␻-frequency trap, the effective value being used for the time evolution is thus scaled by 1/冑² 关see Eq. 共5兲兴, while ␶ is rescaled with respect to the final trap frequency ␻/2. We can see the ‘‘breathing’’ of the conden- sate as it expands and recontracts in the trap.

It is also interesting to look at the effect of the value of the nonlinear parameter on the breathing frequency of the con- densate, as seen in Fig. 3. First, we note that the initial den- sity at the center of the trap is lower for bigger values of␭, which is expected because of the corresponding higher inter-

particle repulsion. Starting from an unperturbed harmonic oscillator (␭⫽0), for which the complete cycle time is ␶

⫽4␲ with recurrences every␶⫽␲, we observe that the os- cillation frequency of the condensate in the trap increases with a greater value of␭.

For the 3D case, we will study the scissor mode关39,40兴 of a trapped condensate. We consider a pancake-shaped con- densate, formed in an anisotropic trap with␻x⫽␻yⰆ␻z, see Fig. 4. The y and z axes of the trap are instantaneously ro- tated, at t⫽0, by an angle␪ around the x axis. The conden- sate then starts to oscillate in the trap, leading to the so- called scissor mode.

FIG. 1. Trapping potentials␻ 共solid line兲 and ␻/2 共dashed line兲 used to simulate the breathing modes of a condensate. The wave function 兩␺(r)兩² of the stationary state for potential ␻ with ␭

⫽100 is also given 共dotted line兲.

FIG. 2. ‘‘Breathing’’ of the condensate after expansion from a trap of frequency␻ to ␻/2. The density profile 兩␺(r)兩²is given as a function of time␶ 共scaled with respect to the final trap frequency

␻/2) for an initial ␭⫽10 ground stationary state.

FIG. 3. ‘‘Breathing’’ of the condensate after expansion from a trap of frequency␻ to ␻/2. The value of the wave function in the center of the trap, 兩␺(r⫽0)兩², is given as a function of time ␶ 共scaled with respect to the final trap frequency␻/2) for ␭ equal to 0 共solid line兲, 10 共dotted line兲, 100 共dashed line兲, and 1000 共dot- dashed line兲.

FIG. 4. Representation of the study a condensate’s scissor mode.

The condensate is initially tilted with respect to the trap’s y and z axes by an angle␪.