Phenomenological damping in trapped atomic Bose-Einstein condensates

5 Pages·0.12 MB·English

by  S. Choi

Phenomenological damping in trapped atomic Bose-Einstein condensates
S. Choi, S. A. Morgan, and K. Burnett
Clarendon Laboratory, Department of Physics, University of Oxford, Parks Road, Oxford OX1 3PU, United Kingdom.
The method of phenomenological damping developed by Pitaevskii for superfluidity near the λ
point is simulated numerically for the case of a dilute, alkali, inhomogeneous Bose-condensed gas
nearabsolutezero. Westudyseveralfeaturesofthismethodindescribingthedampingofexcitations
in a Bose-Einstein condensate. In addition, we show that the method may be employed to obtain
numerically accurate ground states for a variety of trap potentials.
Recentexperiments with dilute, trapped Bose-Einsteincondensates (BEC)haveindicatedthe presenceofdamping
8 ofcollectiveexcitations[1,2]. ThemostcommonlyusedtheoreticaltoolsofarfordescribingBEC,theGross-Pitaevskii
9
equation (GPE) [3], has been shown to give accurate predictions for the frequencies of elementary excitations [4–6],
9
but doesnotcontainanymechanismfor the damping ofthese excitations. Theoriesthatgobeyondthe GPEarenow
1
subjects of intense study [7,8]. One of the theories which incorporates dissipation in a quantum fluid is Pitaevskii’s

n phenomenologicaltheory of superfluidity near the λ point [9], and we have used it to investigatedamping in a dilute,
a
trapped BEC numerically. This phenomenological method of including damping should, as argued by Pitaevskii, be
J
valid when the system is close to thermal equilibrium. We shall begin by deriving the required modification of GPE.
8
We intendtodemonstrateinthis papertwodistinctusesofthe equation: itsuse indescribingdampingofoscillations
2
intrappedatomicBEC,anditsuseasacomputationalmethodtogenerateaccurateeigenstatesofthestandardGPE.

1 The standard GPE describes the motion of the mean field, ψ, which is taken to be the ensemble average of the
v boson field operator Ψˆ. We write the GPE in the form
4
6 ∂ψ ¯h2
0 i¯h =− ∇2rψ+Vtrap(r)ψ+Cn|ψ|2ψ, (1)
∂t 2m
1
0
with C defined by
8 n
9 4π¯h2Na
/ C =NU = , (2)
h n 0
m
p
- where N is the number of particles, a the interatomic s-wave scattering length and m the mass of a single particle.
t
n The mean field ψ is customarily taken to be the condensate wave function, although it is also consistent to define ψ
a to be the condensate ground state plus any coherent excitation.
u
From inspection, the GPE contains no term which can describe damping. The first two terms are simply the
q
energy of a single particle in the trap while the third describes the non-dissipative effect of other particles in the
:
v system. Attempts toimproveonthe GPEinclude takingintoaccountthe correlationeffects ofexcitations[7]andthe
i
X introductionofthe many-bodyT-matrix[8]. Thesetheoriessofarhavenotproducedanyexplicitresultsfordamping
in the inhomogeneous case, although recent related work by Giorgini has shown the presence of Landau and Beliaev
r
a damping [10].
Itis,however,possibletoincludedampingphenomenologically,inamannersimilartothatintroducedbyPitaevskii
[9]. The argument which we present below to obtain our expression for phenomenological damping follows this work
closely. The procedure described is, in fact, a widely used one in the more general problems of dynamical critical
phenomena [11].
Firstofall,wenotethatanydampingprocesseventuallyleadstoanequilibriumstate. ForthecondensateatT =0
such an equilibrium is described by the equation
¯h2
− ∇2ψ+V (r)ψ+C |ψ|2ψ−µψ =0. (3)
r trap n
2m
This expression may be derived from minimizing an appropriate energy functional or equivalently by finding the
eigensolution to the time dependent GPE. Because the process we want to describe is a relaxation process, in our
equation of motion for the condensate ψ,
∂ψ
i¯h =Lˆψ, (4)
∂t
1
the operator Lˆcannot be Hermitian. The anti-Hermitianpart of Lˆis associatedwith the processes by which equilib-
rium is approached. This anti-Hermitian part must vanish at equilibrium, that is, when Eq. (3) is satisfied. For the
cases where the deviation from equilibrium is small, one may therefore write the anti-Hermitian part in the form
¯h2
iΛ ∇2ψ+V (r)ψ+C |ψ|2ψ−µψ , (5)
r trap n
2m
(cid:18) (cid:19)
where Λ is a dimensionless factor which is inversely proportionalto the relaxationtime. The Hermitian part is given
by the fact that for Λ = 0, the theory must reduce to the conventional GPE at T = 0; the final equation including
relaxation to equilibrium is, then, given by
∂ψ ¯h2
i¯h =(1+iΛ) − ∇2+V (r)+C |ψ|2−µ ψ, (6)
r trap n
∂t 2m
(cid:26) (cid:27)
where Λ<0 for damping.
Equation(6)is,infact,arathergeneralequationofmotionwhichdescribesevolutiontowardsequilibrium. Damping
ofelementaryexcitationsmaybe describedintheframeworkofEq. (6)byconsideringourmeanfieldψ tobebroken
into two parts: the condensate ground state ψ and a coherent excitation, δ,
g
ψ =e−iµt(ψ +δ). (7)
g
Evolution in Eq. (6) results in the damping of δ towards zero. The exact value for the damping parameter Λ is
expectedto dependonavarietyoffactors. One would,forinstance,expectit tobe temperaturedependent, since the
thermalcomponent presentinthe trap is clearlygoing to be one ofthe main sourcesofdamping. In the originalcase
ofabsorptionofsoundinsuperfluidHeliumconsideredby Pitaevskii,Λ couldbe expressedinterms ofthe coefficients
ofsecondviscosityusedbyKhalatnikov[12]. Acloselyrelatedresult(whichactuallyderivesthe Landau-Khalatnikov
two-fluidmodelforthe homogeneousBose-condensedgas)hasbeencalculatedbyKirkpatrickandDorfmanusingthe
Chapman-Enskog procedure [13]. They also explicitly derived various transport coefficients such as the viscosity for
their homogeneous case. At present, a complete analysis for the inhomogeneous case is not yet available.
Some rough initial estimate for the damping coefficient Λ may, however, be made by noting that, physically, Λ
represents the rate at which the excited components turn into the condensate. This may be approximated by the
transition probability W+(N) that has been estimated in a recent work by Gardiner et al. [14] using the quantum
kinetic theory. W+(N) givesthe rateatwhich the thermalparticles abovethe condensatebandenter the condensate
due to collisions. Although the description of the process of evaporative cooling was their primary goal, their main
assumption in formulating W+(N) was that the condensate does not readily act back on the thermal component
to change its temperature. This assumption is expected to be valid when one considers quasiequilibrium situations.
Using the reported values from the MIT experiment [2] as an example, and the expression
4m(akT)2 µ µ
W+(N)= e2µ/kT NK N , (8)
π¯h3 kT 1 kT
n (cid:16) (cid:17)o
where K (z) is a modified Bessel function, we obtain Λ∼W+(N)/ω ≈−0.03,for T ≈T /10 and a≈3.45nm. Here,
1 c
ω was taken to be the trap frequency relevant to the oscillatory motion, that is ω = 2π×19Hz. The corresponding
damping time is then approximately 200 ms, which is of the same order of magnitude as the experimental damping
time of 250 ms [2].
Numerical simulation of Eq. (6) is done using a 4th order Runge-Kutta algorithm. Although the simulations in
this paper areperformedusinga onedimensional(1D)Runge-Kuttaintegrator,the generalargumentsthatfollowdo
not depend on dimensionality, and it is straight forward to apply the ideas of the present paper to a 3D simulation.
We simulate damping of elementary excitations by first populating a mode of elementary excitation on a condensate
and then time-evolving this condensate plus excitation using Eq. (6). In practice, such a state is generated by, for
instance,applyingtimedependentperturbationstotheconfiningpotentialinitiallyandthenobservingthesubsequent
behaviourofthecondensateafterturningofftheperturbation. Exactlyhowtheexcitationsdampusingthisformalism
can be analysedquantitatively using the method of quasiparticle projection[15]. This projection method is basedon
the Bogoliubovtransformationof fluctuations aboutthe condensate meanfield: one writes a generalwavefunction ψ
as
ψ =e−iµt (1+b )ψ (r)+ {u (r)b +v∗(r)b∗} (9)
g g i i i i
" #
i
X
2
wheretheindexidenotesquasiparticleenergylevels,andthefunctionsu (r)andv (r)obeyasetofcoupledequations
i i
usually knownas the Bogoliubov-deGennes equations whichdiagonaliseaquadratic Hamiltonian[16]. The sumover
i on the right hand side of Eq. (9) represents the fluctuation about the condensate ground state ψ which is itself
g
weighted by a factor (1+b ). The population of quasiparticles in level i is then given by the quantity b∗b = |b |2
g i i i
where the coefficients b are taken to be the mean values of the quasiparticle annihilation operators. By using the
i
well known orthogonality relations of u (r) and v (r), one may extract the expansion coefficients b and b uniquely
i i g i
at any point in time, and hence deduce the evolution of the quasiparticle population in each mode [15].
In order to simulate the experimental results as closely as possible, we initially populate the breathing mode, the
second lowest excitation, by an arbitrary amount and then plot the evolution of this population, |b |2, over time.
2
Figure 1 shows the populationof mode 2 overtime as well as the populationin the condensate mode |b |2 for various
g
values of Λ. An initial population of 0.25means |b |2 =0.25. The correspondinggroundstate populationis shownin
2
the lower part of the plot. We see that for smaller values of Λ, i.e. when the excited component does not damp out
quickly over time, there is a greater interplay between the ground and the excited state populations due to nonlinear
mixing.
The width ofthe condensateis indeedwhatis actuallymeasuredinthe experiment, andthereforethis is plottedas
afunctionoftimeinFig. 2. ThegeneralshapereportedinRef.[1]oftheformAexp(−Λt)sin(2πνt+φ)+B isclearly
reflected in this figure. The solid curve gives what is expected for the MIT experiment [2] (i.e. Λ ≈ −0.03). As the
simulation was done in 1D, however, we do not intend to claim any quantitative analogy. Also, in real experiments,
several quasiparticle modes, rather than just the breathing mode are simultaneously excited. It is clear, however,
that this descriptionofdamping shouldbe satisfactoryincaseswhere anexactformalismis not needed: forinstance,
when Λ may be treated as a parameter to be obtained from an experiment.
One of the ways to characterise different descriptions of damping is to see its dependence on energy. For instance,
in the homogeneous gas, it is well known that the damping rate increases as k2 when the relaxation process is due
to thermal conduction and viscosity [11]. We note that the Pitaevskii prescription gives a damping rate which has
a linear dependence on the mode energy, with the high energy components being damped out more rapidly. This
aspect may be seensimply by inspection ofthe equation ofmotion Eq.(6), and it can alsobe confirmed by explicitly
calculating the damping rate of the quasiparticle populations, |b |2, for various values of i by linearisation of Eq. (6).
i
We nowlookatthe secondpossibleuse ofEq.(6). We notethatfinding areliablesolutiontothe time-independent
GPE is an essential first step in many numerical simulations involving a trapped BEC. The Pitaevskii damping
description may, in fact, be used to this end as an effective numerical tool to find the eigenstate. This is possible
because the equilibrium state Eq. (3) is itself a solution to Eq. (6). To find the eigenstate, one starts by initially
choosing an arbitrary function (preferably close to the desired solution for faster convergence, e.g. for a condensate
in a trap, a Gaussian or a Thomas-Fermi solution) and run the routine with pre-determined values of µ and Λ. It
should be noted here that, since the equation is non-hermitian, the evolution does not conserve the norm of ψ. It is
thereforenecessaryto proceedby firstrunning the routine with anestimate of C forthe chosenµ until convergence,
n
and then later adjusting C from the renormalisation of ψ,
n
|ψ(r)|2dr=1. (10)
Z
The correct C for the chosen µ is then given by the initial estimate divided by the final normalisation constant.
n
We found that the time it takes to find a solution in this way is very much shorter than other methods such as the
adiabaticturn onof the non-linearityC . As aroughindication,evolutiontime oft=10/ω wasrequiredforgrowing
n
an eigenstate with C = 100 and Λ = −2 when using the damping method while it took us t = 500/ω using the
n
adiabatic method; in real time on a Sun SPARC station 10, these running times corresponded to roughly 10 minutes
and 3 hours respectively. The damping method is able to find a solution for as large a C as desired, provided there
n
are a sufficient number of grid points to accommodate the larger condensate. An exact eigensolution to GPE has a
flat phase profile; the damping simulation with Λ = −2, C = 50 running for t = 15/ω typically gave a maximum
n
change in phase across the condensate of the order of 10−7 radians, indicating the high accuracy of this method.
Other advantages of this approach include its flexibility in handling traps of any arbitrary shape, and the fact
that it is not necessary to expand the wave function in a set of complete, orthogonal basis states which can be a
computationally expensive operation.
The value of Λ cannot, however, be increased indefinitely for faster convergence, as this results in numerical
instability. ItwasfoundthatwegetnumericalinstabilitywhenΛt >0.03wheret isthe sizeofthe incremental
step step
time step in numericalintegration. The onsetof numericalinstability is independent ofthe total lengthof simulation
andthe sizeofthe non-linearityconstantC . Onemayunderstandthis fromthe factthatthe dampingper time step
n
is roughly given by e−Λtstep so that e−0.03 ≈ 0.97 implies a reduction in the amplitude by about 3 % per numerical
time step. Since the typical step size used is of the order of 10−3 harmonic oscillator units, this, in fact, represents a
complete decay on a time scale of order a period.
3
We also point out that the vortex eigenstate in a trap may be obtained numerically by using this damping pre-
scription; one needs only to start with the first excited harmonic oscillator state which has an odd parity. One of the
corollariesofthedampingmethodisthatalthoughparityisconservedthenumberofzerocrossingsisnot;ifwestarted
with, say,the3rdexcitedharmonicoscillatorstate,the resultingsolutionisstillthe statewithonezerocrossing. The
phase profile of such a condensate displays a clear π phase jump at the origin, while being flat everywhere else.
Insummary,wefoundthatthephenomenologicaldampingdescriptiondoesindeedgiveresultswhicharecompatible
with those of the experiments. In a more general context, it provides an efficient numerical tool for finding the
eigenstate of the time independent GPE.
ACKNOWLEDGMENTS
SC would like to acknowlegethe support ofUK CVCP andthe Wingate Foundation; SM and KB acknowledgethe
support from the UK EPSRC. KB also acknowledges support from the European Union under the TMR Network
Programme.
[1] D.S. Jin, M. R.Matthews, J. R. Enscher, C. E. Wieman, and E. A. Cornell, Phys. Rev.Lett. 78, 764 (1997)
[2] M.-O.Mewes, M.R.Andrews,N.J.vanDruten,D.M.Kurn,D.S.Durfee,C.G.Townsend,andW.KetterlePhys.Rev.
Lett.77, 988 (1996)
[3] V. L. Ginzburg and L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 34, 1240 (1958) [Sov. Phys. JETP 7, 858 (1958)]; E. P. Gross,
J. Math. Phys. 4, 195 (1963)
[4] M. Edwards, P.A. Ruprecht,K. Burnett,R. J. Dodd, and C. W. Clark Phys.Rev.Lett 77, 1671 (1996)
[5] D.S. Jin, J. R.Enscher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Phys. Rev.Lett. 77, 420 (1996)
[6] S.Stringari, Phys.Rev. Lett 77, 2360 (1996)
[7] A.Griffin Phys. Rev.B 53, 9341 (1996)
[8] N.P. Proukakis, K. Burnett,and H.T. C. Stoof, To appear in Phys. Rev.A (1998)
[9] L. P.Pitaevskii, Sov. Phys.J.E.T.P. 35, 282 (1959)
[10] S.Giorgini, To appear in Phys.Rev.A (1998), e-print cond-mat/9709259
[11] E.M.LifshitzandL.P.Pitaevskii,PhysicalKinetics,LandauandLifshitzCourseofTheoreticalPhysicsVol.10(Pergamon
Press, Oxford 1981)
[12] I.M. Khalatnikov J. Exptl. Theoret. Phys. (U.S.S.R.) 23, 8 (1952)
[13] T. R. Kirkpatrick and J. R. Dorfman, Phys. Rev.A 28, 2576 (1983)
[14] C. W. Gardiner, P. Zoller, R. J. Ballagh, and M. J. Davis, Phys. Rev.Lett. 79, 1793 (1997)
[15] S.A. Morgan, S. Choi, K. Burnett, and M. Edwards submitted toPhys. Rev.A
[16] A.L. Fetter, Ann.Phys.70, 67 (1972)
4
Evolution of quasiparticle and condensate population
0.3
0.25 quasiparticle
0.2
0.15
0.1
population 0.05
0
−0.05
−0.1
−0.15 condensate
−0.2
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
time in h.o. periods
FIG.1. Populationofthebreathingmode,|b2|2,overtimeforvariousvaluesofΛ. Solid,dashed,dot-dashed,anddottedlinesindicate
Λof−0.03,−0.1,−0.25,and−0.5respectively. Thelowerhalfoftheplotshowsthecorrespondingvariationin|bg|2.
Evolution of condensate width over time
2
1.8
1.6
width of condensate in h.o. units11..24
1
0.8
0.6
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
time in h.o. periods
FIG. 2. Variation of the width of the trapped condensate over time. The width is measured by the standard deviation of x, the
position,weightedbythecorrespondingcondensatedensity. Solid,dashed,dot-dashedanddottedlinesindicateΛof−0.03,−0.1,−0.25,
and−0.5respectively,asinthepreviousfigure.
5