18 October 1999
Physics Letters A 261 Ž1999. 337–344 www.elsevier.nlrlocaterphysleta
Excitation frequencies and stability regions for two-species Bose–Einstein condensate in a triaxial magnetic trap A. Montina, R. Mannella, E. Arimondo Dipartimento di Fisica and Istituto Nazionale Fisica della Materia, UniÕersita` di Pisa, Via Buonarroti 1, 56100 Pisa, Italy Received 2 August 1999; accepted 9 August 1999 Communicated by V.M. Agranovich
Abstract We have solved the coupled Gross–Pitaevskii equations for a two-species condensate confined within a fully anisotropic magnetic trap. Both variational and numerical methods are used to solve the Gross–Pitaevskii equations, and the original numerical approach is discussed in more length. The calculation of the two-species excitation spectrum shows that the mode frequencies of the individual condensates are modified in a two-species condensate interatomic interaction, with new excitation modes present in the mixed system. The spatial distributions of the two condensates show deformations produced by the interatomic interactions, with either compression of the atoms in the central region of the trap, or with a decrease of the atomic density in that region. q 1999 Published by Elsevier Science B.V. All rights reserved. PACS: 03.75.Fi; 05.30.Jp; 34.50.Rk; 33.70.C
1. Introduction In the framework of Bose–Einstein condensation, the observation of a condensate mixture composed of two spin states of 87 Rb w1x has prompted a significant interest into the two-species Bose–Einstein condensates ŽBEC’s.. In these systems two different atomic species are simultaneously confined in a magnetic trap and cooled to very low temperatures, where most atoms are described by the same macroscopic wavefunction. In the condensation phase, the atomic properties are determined by the confining potential of the magnetic trap, but also by the collisions within each atomic species. For the two-species condensate the interatomic collisions play an addi-
tional role. The presence of an additional term in the condensate interaction energy produces an extra freedom in the construction of the Bose–Einstein condensate. Several important properties of the twospecies Bose condensate have been previously examined w2–10x, the most important one being the spatial symmetry breaking associated to the macroscopic wavefunction of the condensate. Thus in a symmetric magnetic potential the condensate wavefunction does not match to the trap symmetry w4x. Also, the phenomenon of density compression may take place, with one atomic species squeezed out and forced to form a shell around the trap center which holds the second atomic species, receiving extra confinement from the first atomic species w6x.
0375-9601r99r$ - see front matter q 1999 Published by Elsevier Science B.V. All rights reserved. PII: S 0 3 7 5 - 9 6 0 1 Ž 9 9 . 0 0 5 7 6 - 9
338
A. Montina et al.r Physics Letters A 261 (1999) 337–344
In this work we studied a physical system composed by two Bose–Einstein condensates, rubidium and cesium, confined in a triaxial fully anisotropic magnetic trap. Throughout the paper we used the well known rubidium-rubidium scattering length, and the current best estimates for the cesium-cesium and rubidium-cesium scattering lengths. The estimated large negative scattering length of cesium implies that a condensate of cesium atoms is stable only for a very small number of bosons. Our aim is to investigate whether, using the density compression provided by a rubidium condensate obtained using a large number of bosons, it is possible to realize a stable condensate of cesium atoms with a sizable number of atoms, i.e., larger than that achievable in presence of single species. This should give some indications whether, as a final product of the sympathetic cooling between two different alkali species, as opposed to two isotopes of the same species, two species condensates could be created even in very unfavorable cases. The starting point of our two-species BEC investigation is the calculation of the two-species excitation spectrum. We show that the mode frequencies of the individual condensates are modified by the interatomic interaction and that new excitation modes are present in a two-specie condensate. To deal with a condensate containing a small number of atoms, like cesium with a negative scattering length, the Thomas–Fermi approximation cannot be used, so we used both a variational approach based on Gaussian trial functions and a direct numerical integration of the tridimensional time dependent GPE. The numerical solution allowed us to derive the deviation in the spatial distribution from a Gaussian shape for both atomic species, with either a compression or a depression at the center of the magnetic trap. Our analysis seems to support the conjecture that, even in the presence of very few cesium atoms, the perturbation induced on the rubidium atoms by the presence of the second boson should lead to observable effects, both in the spatial distribution of the rubidium cloud and in the collective frequencies. The organization of the manuscript is based on the initial presentation of the variational and numerical methods for the GPE solution. Then, our results for the spatial distribution of the condensates and the collective frequencies for the condensate oscillations
are examined as a function of the number of atoms. Conclusions and some final considerations are drawn in the closing section.
2. GPE and variational method In mean field approximation the Bose condensate with two species is described by the macroscopic wave functions c 1 and c 2 , that evolve in time according to the self-consistent Gross–Pitaevskii equations ŽGPE’s. E i"
Et
c 1Ž r ,t . " 2=
ž
2
q V 1 Ž r . q g 11 < c 1Ž r ,t . < 2 q g 12 < c 2 Ž r ,t . < 2 c 1Ž r ,t . ,
Ž 1a .
q V 2 Ž r . q g 22 < c 2 Ž r ,t . < 2 q g 12 < c 1 Ž r ,t . < 2 c 2 Ž r ,t . ,
Ž 1b .
/
s y 2 m1
E i"
Et
c 2 Ž r ,t . " 2=
ž
2
s y 2m2
/
where m1 and m 2 are the masses of the two atoms. The coupling constants g 11 , g 22 g 12 are given by g i i s 4 p " 2 a i irm i and g 12 s 2 p " 2 a12 Ž m 1 q m 2 .rm1 m 2 , being a i i e a12 respectively the s-wave scattering lengths between atoms of the same species and of different species. Eqs. Ž1. are formally derived minimizing the action S defined by S s L d 3 r dt
Ž 2.
H
associated to the Lagrangian density L s i" c 1) y
Ec 1 Et
g 11 2
"2 y 2 m1
=c 1)=c 1 y V1 Ž r . < c 1 < 2
< c 1 < 4 q i" c 2)
y V2 Ž r . < c 2 < 2 y
g 22 2
Ec 2 Et
"2 y 2 m2
=c 2)=c 2
< c 2 < 4 y g 12 < c 1 < 2 < c 2 < 2 .
Ž 3. Like in the BEC experiments in magnetic traps for single atomic species, we have assumed that the two condensates are spatially confined inside a magnetic trap with the following triaxial harmonic potential Ža triaxial TOP trap as in the experiments of Ref. w11x. 2 V1 Ž r . s 12 m1 v 01 Ž l12 x 2 q l22 y 2 q l23 z 2 . ,
V2 Ž r . s
1 2
2 m 2 v 02
Ž
l12 x 2 q l22
y
2
q l23 z 2
..
Ž 4a . Ž 4b .
A. Montina et al.r Physics Letters A 261 (1999) 337–344
Here the adimensional constants l i , scaling the spring constants along the three axes, were assumed to be 1, '2 , and 2 respectively. We supposed that the gravity sag of the two condensate are exactly compensated in order to produce the two harmonic potentials centered the same point in space w1x. To solve Eqs. Ž1. we used two different approaches, detailed below: a variational method based on trial Gaussian wave functions, and a direct numerical solution of the 3 q 1 GPE equation. As trial functions in the variational analysis we have used
c 1 Ž r ,t . s N11r2
ž
1r4
1
p 3s 112 Ž t . s 212 Ž t . s 312 Ž t .
/
Ł
is1,2,3
2
½
=exp y
Ž x i y hi1 . q i a i1 x i q i bi1 Ž t . x i2 , 2 si12 Ž t . Ž 5a .
5
339
where U Ž s ,h .
(
s4
a12 " 2 N2 Ž m1 q m 2 .
1
p m1 m 2
(Ł Ž k
s k12 q s k22 .
Žh k1 y h k 2 .2
=ey Ý s k1 2q s k 2 2 .
Ž 7.
k
Exchanging the indexes 1 and 2 in Eq. Ž6. we obtain similar equations for the parameters of the second condensate species. Note that the quantities a i j and bi j do not enter Eq. Ž6.. The equilibrium configurations for the system are found imposing that the time derivatives of s and h are equal to zero. Linearizing the equations of motion around the equilibrium position and diagonalizing them we have obtained the collective oscillation frequencies for the two-species condensate.
c 2 Ž r ,t . s N21r2
ž
1r4
1
p 3s 122 Ž t . s 222 Ž t . s 322 Ž t .
/
3. Numerical solution
Ł
is1,2,3
2
½
=exp y
Ž x i y hi2 . q i a i2 x i q i bi2 Ž t . x i2 , 2 si22 Ž t . Ž 5b .
5
where si j and hi j are respectively the widths and positions of the two Gaussians in the three spatial directions x i , while a i j and bi j are related respectively to the translation and expansion rates of the distributions. Minimizing the action S in Eq. Ž2. we obtain the following Newton-like equations for the quantities si1 and hi1:
Esi1
k
Et
ž
c Ž x ,t . s y
"2 E 2 2m E x2
/
q V Ž x . c Ž x ,t .
' H Ž x , c ,t . c Ž x ,t . .
Ž 8.
"2
E
qŁ
E i"
This equation can be integrated by a sort of Crank– Nicholson scheme w13x
2 s¨i1 Ž t . q l2i v 01 si1
sy
For the numerical integration of the time dependent GPE equation we have used a modified split operator technique, adapted to the integration of a Schrodinger equation, this algorithm already briefly ¨ accounted in w12x, but here opportunely extended to the two species case. To see how the algorithm works, we start writing a 1 q 1 Schrodinger equation ¨
(
Ý 2 m2 s 2 Ž t . 1
k
2 a11 " 2 N1
p
m12 s k1
E 2 h¨ i1 q l2i v 01 hi1 s y
ž
k1
q
N2 m1
N2
Ehi1 m1
U Ž s ,h .
U Ž s ,h . ,
Ž 6a . Ž 6b .
1q
iHd 2"
/
ž
c Ž x ,t q d . s 1 y
iHd 2"
/
c Ž x ,t . ,
Ž 9.
where d is the integration time step. Discretizing this equation on a lattice, due to the second order spatial derivative on the left hand side, it is clear that each integration time step requires the solution of a tridiagonal linear system. To keep the following
A. Montina et al.r Physics Letters A 261 (1999) 337–344
340
algebra simple, we rewrite Eq. Ž9. in the equivalent form
ž
c Ž x ,t q d . s 1 q
iHd 2"
y1
/ ž
1y
iHd 2"
/
c Ž x ,t . .
Ž 10 . Now, let us go back to the original 3 q 1 GPE. Clearly, we cannot apply the same algorithm here. The discretization on a lattice of the Žcorresponding. left hand side would lead to a sparse linear system not in tridiagonal form. In particular, having discretized each spatial direction with N bins, we would have to solve a N 3 = N 3 complex linear system at each integration time step, although only 7N y 6 elements on the left hand side are different from zero. The idea is then to split the 3 q 1 GPE Hamiltonian, and integrate three equivalent 1 q 1 problems. We start writing the 3 q 1 GPE in the following form Žwe consider here the case of one boson, the extension to more than one boson is simple but cumbersome.:
E i"
Et
c Ž r ,t .
s H x Ž r ,t . q H y Ž r ,t . q Hz Ž r ,t . c Ž r ,t . ,
Ž 11 . where "2 E 2 Hi Ž r ,t . ' y
2 m E x i2
q V Ž x i . q 13 g < c Ž r ,t . < 2 .
Ž 12 . The splitting is carried out so that the commutators are exact up to the order d 2 included. Eq. Ž11. was integrated using the following scheme, where d is the integration time step and A i Ž t . ' i d Hi Ž r,t .r" .:
c Ž r ,t q d . s 1 q A y Ž t . r2
y1
1 y A x Ž t . r2
= 1 q A z Ž t . r2
y1
1 y A z Ž t . r2
= 1 q A x Ž t . r2
y1
1 y A y Ž t . r2
=c Ž r ,t . .
Ž 13 .
To be more explicit, at each integration time step we apply first 1 y A y Ž t . r2 to c ; then we solve the
tridiagonal system we obtain with 1 q A x Ž t . r2 on the left hand side and get an intermediate c , we apply 1 y A z Ž t . r2 to this intermediate c and so on. There is obviously a problem with the nonlinear term g < c Ž r,t .< 2 , because we should really use a c somehow averaged over the time step d , not a c evaluated at the beginning of the time step. To circumvent this problem, we used a sort of predictor corrector step. Each integration step is really done in two times: going from the time t to the time t q d , the first time we use c Ž r,t . in the nonlinear term, obtaining a ‘‘predicted’’ c Ž r,t q d .; we then repeat the integration step, starting again from c Ž r,t ., but using 12 Ž c Ž r ,t . q c Ž r ,t q d . . in the nonlinear term. To find the steady state we evolved the condensate wavefunctions in complex time Žsee also w14x.. However, contrary to what is done in w14x, we did not use the chemical potential. In our numerical solution of the GPE, we fixed the particle number for the two species and we renormalized the wavefunctions at each integration time step. It is clear that the steady state of the GPE Ž c 0 . is also a steady state for our algorithm: after an integration step in complex time, c 0 is modified by a numerical factor independent of the spatial coordinates; after applying the renormalization we have again c 0 apart from an overall Žirrelevant. phase factor. The oscillation frequencies of the condensates were found applying a small perturbation to the steady state, evolving the condensate in real time, and looking at the frequency spectrum of the widths and positions of the two condensates.
4. Results and discussion In the numerical simulations we have considered atoms 133 Cs and 87 Rb, trapped respectively in the < F s 3, MF s y3: and < F s 1, MF s y1: hyperfine sublevels. For the rubidium-rubidium scattering length we used aRb – Rb s 109.1 a.u. for the < F s 1, MF s y1: hyperfine sublevel. For aCs – Cs we have studied both cases of repulsive positive and negative attractive values, assuming that, due to the presence of proper Feshbach resonances, the sign of the scat-
A. Montina et al.r Physics Letters A 261 (1999) 337–344
Fig. 1. Widths of the condensate Gaussian distributions computed with the variational approach Žsolid line. and numerically Žsymbols., as function of the number of bosons of the two condensates Žassumed equal.. The frequencies for the cesium and rubidium magnetic traps are respectively l i v Rb and l i v Cs , with l i s 1,'2 ,2 and v Rb s 2p =565 Hz and v Cs s 2p =559 Hz.
tering length could be modified also for the magnetically trapped < F s 3, MF s y3: hyperfine sublevel. We have used for the scattering length the absolute value of 400 a.u., in agreement with different recent estimates w15x, with the appropriate sign depending on the case studied. Finally, for the rubidium-cesium collisions we have used the estimate aRb – Cs s 150 a.u. 1. This value is smaller than the critical value c < aRb < < < above which in the – Cs s a Rb – Rb a Cs – Cs Thomas–Fermi regime the two-species condensate cannot exist w17,5,9x. We have initially supposed for cesium a positive scattering length in order to test the validity of our approaches. We took the boson number N to be the same for the rubidium and cesium condensates. Fig. 1 shows the widths si j of the Gaussian distributions obtained as a function of the boson number N using the variational and numerical methods, assuming for cesium scattering length the 400 a.u. value and with the magnetic trap parameters given in the figure caption. For a small number of bosons the kinetic energy is larger than the interaction energy, and the rubidium spatial distribution is wider than that of
341
cesium because of its smaller mass. At large N ’s, in the Thomas–Fermi regime where the kinetic energy is smaller than the interaction energy, the situation of the spatial distribution is inverted: the cesium larger scattering length yields a wider Gaussian shape. Increasing N we note that the agreement between the two methods remains good for the cesium data, while a disagreement appears for the rubidium data. The reason of this disagreement is clear from the results reported in Fig. 2, where the probability distributions < c Cs Ž x,0,0.< 2 for cesium, computed numerically and via the Gaussian ansatz, are plotted. The numerical distribution appears, close to the center of the trap, flatter than the Gaussian distribution calculated with the variational method. Therefore this last approach over-estimates the repulsive effect of the cesium atoms on the rubidium ones, and consequently it predicts a larger width for the rubidium cloud. The results of Fig. 2 evidence another important result of our numerical analysis: the interaction between the two alkali species may produce large deviations from the Gaussian distribution. Fig. 3 shows our results for the expected oscillation frequencies. The lines are the prediction from the variational analysis, the symbols are the results of the numerical integration of the GPE. The figure shows two different sets of curves Žsee Eq. Ž6..: the dashed Ždotted. curves are relative to the motion of
(
1 The rubidium-cesium scattering length has been derived on the basis of published spectroscopy data for the ground state of the RbCs molecule, see w16x.
Fig. 2. < c Cs Ž x,0,0.< distribution of the cesium condensate Žsupposing a positive scattering length. in the presence of a rubidium condensate with equal number of atoms, at Ns 1000. The dotted line shows the results from the variational approach and the solid line those evaluated numerically evaluated.
342
A. Montina et al.r Physics Letters A 261 (1999) 337–344
Fig. 3. Collective frequencies evaluated with the variational method Žsolid line. and numerically Žsymbols.. Parameters as in Fig. 1.
the s ’s Žh ’s.. The short solid horizontal lines visible on the right and on the left of the figure are notable limits, which will be described below. The magnetic trap oscillation modes of the center of gravity are for rubidium at vrv x ,Rb s 1,1.414,2 and for cesium at vrv x ,Rb s 0.989,1.399,1.979 2 . For very small N ’s, we expect the frequencies for the breathing motion Ž s ’s. of the condensate Gaussian width to be exactly twice the above ones, hence slightly different for the two species: however, on the vertical scale of the figure the difference between rubidium and cesium cannot be distinguished. Still for small N ’s, we expect the frequencies for the motion of the center of mass of the condensates to be equal to the trap frequencies: this is clearly the case. Increasing the number N of boson atoms, the nonlinear character of the GPE manifests itself: it is found that the breathing collective motion frequencies should go to the hydrodynamic limit Žshown on the right hand side of the figure.. In this limit, the frequencies of a single condensate no longer depend on the scattering length, but only on the trap frequency w18x; because of the chosen magnetic trap configuration Žequal frequencies for the two species.
2 Making use of the different hyperfine levels for the two alkalis it is possible to realize a magnetic trap configuration with the trap frequencies nearly equal for the two atoms.
we find that in our case we should still have some oscillations at frequencies equal to the hydrodynamic limit. On the other hand, it is possible to have a motion where the breathing modes of the two species are ‘‘out of phase’’, which should lead to oscillations frequencies which depend on the scattering lenghts. We see that our variational analysis does reproduce this intuitive pictures, with three breathing mode frequencies going to the hydrodynamic limit. We also note that these frequencies are in good agreement with the frequencies found in the numerical integration of the GPE. There is not agreement between three breathing mode frequencies Žnotably, the ones which should depend on the scattering lenghts. and the result of the numerical integration of the GPE. To assess the origin of this discrepancy, we extended the work of w8x to the triaxial case: it turns out that the numerical solution of the GPE is correct, and that the three frequencies which are not reproduced by the variational method correspond to higher order modes which, as N is increased, are pushed towards lower frequencies. As for the motion of the center of mass of the condensate in the limit of large N ’s, we expect to find the unperturbed trap frequencies Žcorresponding to an in phase oscillation of the two condensate, marked by the solid line on the right of the figure., and some additional frequencies which will depend on the scattering lengths, corresponding to an out of phase motion. This is indeed the case, as confirmed by the numerical solution of the GPE. We turn now to the case of a negative scattering length for cesium. If rubidium is not present, there exists a critical number Nc for the bosons number above which cesium collapses because of the attractive term of auto interaction within the GPE w18x. That critical number can be increased reducing the spring constants of the harmonic magnetic confining potential. Such dependence suggests that condensate stability with a number of bosons greater than Nc can be realized adding a second condensate that interacts repulsively with cesium. This configuration is realized for a negative cesium-cesium scattering length and a positive rubidium-cesium scattering length. Using the variational method we have derived the steady state of the GPE’s by varying NCs and NRb . For NRb different from zero we have found that, increasing NCs , a steady state exists until a critical
A. Montina et al.r Physics Letters A 261 (1999) 337–344
Fig. 4. Stability region for a cesium-rubidium mixture. To the left of the solid line Žobtained from the variational approach. there is the region where the cesium condensate is stable. The circles Ždaggers. mark the values of the cesium and rubidium atom numbers yielding an unstable Žstable. cesium condensate: the condition of stability is derived numerically.
value of the cesium atomic number is reached where the lowest oscillation frequency of the two-species condensate goes to zero. This particular frequency is related to the vibration of the center of mass along one of the axes, and its change in sign indicates that the state has become unstable by broken symmetry of the overall GPE solution 3. In these conditions the two condensates separate spatially moving into two distinct regions, so that they are no longer at the center of the trap. The broken symmetry value coincides with the number of cesium atoms for which the cesium condensate collapses: as soon as the twospecies condensate separates in space, the cesium condensate is no longer sustained against collapse by the other condensate species. Whence for a broken symmetry with a condensate cloud having a number of bosons larger than the one-species critical value, a condensate collapse occurs. Fig. 4 shows the stability region of the two-species condensate on the NCs –NRb plane. The solid line divides the stable region, at lower atom number values, from the unstable region, at larger atom number values, as computed using the variational
3 For a fully anisotropic trap, in the presence of a broken symmetry, the wavefunction invariance under reflection on one of the planes Ž is 0. with Ž is x, y, z . no longer applies.
343
method. We find that, increasing NCs , the critical number of Cs atoms for which we observe the broken symmetry point is larger than the critical Nc value for a single species. Points computed using the numerical integration on the lattice are shown by symbols: values of Ž NCs , NRb . supporting a stable or unstable mixture are shown as diamonds or daggers, respectively. The transition line between stable and unstable mixtures lays between the diamonds and daggers of the figure. It is clear from the figure that the variational method over-estimates the region of stability. The reason of this discrepancy between the two methods is similar to that reported in the discussion of the results of Fig. 1, and evidenced by the spatial distribution of Fig. 2. Fig. 5 shows the probability distribution for the rubidium atoms along the x axis, < c Rb Ž x,0,0.< 2 , obtained by the numerical method in the case of negative aCs ,Cs : at larger NCs , we find that at the trap center instead of a maximum the rubidium distribution has a depression caused by the presence of cesium. This deformation of the rubidium cloud produces an additional squeezing of the cesium distribution towards the trap center. The numerical analysis allows us to conclude that this change in the rubidium curvature at the trap center has two important effects: first, one of the hydrodynamics frequencies connected to the motion of the s Gaussian width goes to zero; second, the collapse of
Fig. 5. Probability distribution < c Rb Ž x,0,0.< of the rubidium condensate in the presence of a cesium condensate with negative scattering length, for the following values of atom numbers, NRb s 6000 and NCs s 8. The line shows results from the numerical approach.
344
A. Montina et al.r Physics Letters A 261 (1999) 337–344
the cesium condensate does not happen via broken symmetry, as for the case of the variational analysis
scattering length. Up to now double species condensates have not yet produced, so that the phenomena here predicted could not be directly tested.
5. Conclusions References We have performed numerical analyses for a double Bose–Einstein condensate confined within a triaxial magnetic trap, using both a variational approach and the numerical solution of the two-species GPE. We have examined several characteristic BEC quantities to be tested in experiments on two-species. The most important features derived here are associated to the spatial separation of the two atomic species, produced the positive interatomic scattering length. For instance at large number of rubidium atoms the interatomic repulsion is so large that the potential experienced by the cesium atoms is no more harmonic. Viceversa the cesium atoms, even at small number, produce a large deformation in the rubidium cloud. It may be imagined that the presence of the small condensate cesium could be detected not directly on the cesium absorption, but on the spatial deformation of the rubidium cloud. In reality for the values of cesium-cesium and scattering lengths we have introduced in our analysis, the modifications in the rubidium cloud are very small and not easily monitored on the integrated absorption profile. However for other scattering lengths, and whence for different atomic species, larger deformations may be produced and detected. The excitation frequencies of a double BEC may be used, as already for the single species BEC, to test the precise knowledge of the interatomic interactions. The collective atomic frequencies depend on the number of atoms present in the two species. For a negative cesium scattering length, where the condensate will contain a critical maximum atomic number, the modification in that critical number could be monitored through the modifications in the collective frequencies going to zero near the collapse of the condensate having negative
w1x C.J. Myatt, E.A. Burt, R.W. Ghrist, E.A. Cornell, C.E. Wieman, Phys. Rev. Lett. 78 Ž1997. 586. w2x Tin-Lun Ho, V.B. Shenoy, Phys. Rev. Lett. 77 Ž1996. 3276. w3x E. Goldstein, P. Meystre, Phys. Rev. A 55 Ž1997. 2935. w4x B.D. Esry, Chris H. Greene, James P. Burke Jr. and John L. Bohn Phys. Rev. Lett. 78 Ž1997. 3594; B.D. Esry, C.H. Greene, Phys. Rev. A 57 Ž1998. 1265; B.D. Esry, ibidem A58 Ž1998. R3399; J.P. Burke Jr., J.L. Bohn, B.D. Esry, C.H. Green, Phys. Rev. Lett. 80 Ž1998. 2097, and Phys. Rev. A 59 Ž1999. 1457. w5x C.K. Law, H. Pu, N.P. Bigelowm, J. Eberly, Phys. Rev. Lett. 79 Ž1997. 3105. w6x H. Pu, N.P. Bigelow, Phys. Rev. Lett. 80 Ž1998. 1130, 1134. ¨ w7x P. Ohberg, S. Stenholm, Phys. Rev. A 57 Ž1998. 1272. w8x D. Gordon, C.M. Savage, Phys. Rev. A 58 Ž1998. 1440. w9x P. Ao, S.T. Chui, Phys. Rev. A 58 Ž1999. 4836. w10x A. Sinatra, P.O. Fedichev, Y. Castin, J. Dalibard, G.V. Shlyapnikov, Phys. Rev. Lett. 82 Ž1999. 251. w11x L. Deng. E.W. Hagley, J. Wen, M. Trippenach, Y. Band, P.S. Jiulienne, J.E. Simsairan, K. Helmerson, S.L. Rolston, W.D. Phillips, Nature 398 Ž1999. 218; E.W. Hagley, L. Deng, M Kozuma, J. Wen, K. Helmerson, S.L. Rolston, W.D. Phillips, Science 283 Ž1999. 1706; Yu. B. Ovchinnikov, J.H. Muller, M.R. Doery, E.J.D. Vredenbregt, K. ¨ Helmerson, S.L. Rolston, W.D. Phillips, Phys. Rev. Lett. 83 Ž1999. 284. Ž1983.; 2082 Ž1991.. M. Lewenstein. w12x E. Cerboneschi, R. Mannella, E. Arimondo, L. Salasnich, Phys. Lett. A 249 Ž1998. 495. w13x W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes. The Art of Scientific Computing, Cambridge University Press ŽCambridge, 1989.. w14x S. Choi, S.A. Morgan, K. Burnett, Phys. Rev. A 57 Ž1998. 4057. w15x V. Vuletic, ´ A.J. Kerman, C. Chin, S. Chu, Phys. Rev. Lett. 82 Ž1999. 1406, and references therein. w16x A. Bambini, private communication. w17x E.D. Sigga, A.E. Ruckenstein, Phys. Rev. Lett. 44 Ž1980. 1423; A.F. Andreev, E.P. Bashkin, Zh. Eksp. Teor. Fiz. 69 Ž1957. 319 wSov. Phys. JETP 7 Ž1958. 698x. w18x F. Dalfovo, S. Giorgini, L.P. Pitaevskii, S. Stringari, Rev. Mod. Phys. 71 Ž1999. 463.