Ultra-Fast Converging Path-Integral Approach for Rotating Ideal Bose-Einstein Condensates

1.6 MB·English

by  Antun Balaz

Ultra-Fast Converging Path-Integral Approach for Rotating
Ideal Bose-Einstein Condensates
Antun Balaˇz∗,a,1, Ivana Vidanovi´ca,1, Aleksandar Bogojevi´ca,1, Axel Pelsterb,c
aScientific Computing Laboratory, Institute of Physics Belgrade, Pregrevica 118, 11080 Belgrade, Serbia
0 bFachbereich Physik, Universita¨t Duisburg-Essen, Lotharstraße 1, 47048 Duisburg, Germany
1 cUniversita¨tPotsdam, Campus Golm, Karl-Liebknecht-Strasse 24/25, 14476 Potsdam-Golm, Germany
0
2

n
a
J Abstract

1
1 A recently developedefficient recursiveapproachfor analytically calculating the short-time evolutionof the
one-particle propagator to extremely high orders is applied here for numerically studying the thermody-

] namical and dynamical properties of a rotating ideal Bose gas of 87Rb atoms in an anharmonic trap. At
s
a first, the one-particle energy spectrum of the system is obtained by diagonalizing the discretized short-time
g
propagator. Using this, many-boson properties such as the condensation temperature, the ground-state
-
t occupancy, density profiles, and time-of-flight absorption pictures are calculated for varying rotation fre-
n
quencies. The obtainedresults improve previoussemiclassicalcalculations,in particularfor smaller particle
a
u numbers. Furthermore, we find that typical time scales for a free expansion are increased by an order of
q magnitude for the delicate regime of both critical and overcriticalrotation.
.
t
a Key words: Bose-Einstein condensation, Effective action, Exact diagonalization,Time-of-flight
m PACS: 03.75.Hh, 31.15.xk,51.30.+i
-
d
n 1. Introduction to expand perpendicular to the rotation axis. For
o
an overcritical rotation, this expansion would even
c Bose-Einstein condensation (BEC) represents a
[ beacceleratedbythepresenceofaresidualcentrifu-
macroscopicquantum phenomenonof recent broad
galforce. Inordertoreachexperimentallythisdel-
1 research interest [1, 2, 3, 4]. Since its first experi-
icate regime of criticalor overcriticalrotation, Fet-
v mental realization in 1995, it has been extensively
3 tersuggestedinRef.[6]toaddanadditionalquartic
studied experimentally, analytically, and numeri-
6 termtotheharmonictrappotential. UsingaGaus-
cally. Thetwomainresearchdirectionsareweakly-
4 sianlaserbeampropagatinginthe z-direction,this
1 interacting dilute gases in magneto-optical traps
has been realized experimentally in Paris by Dal-
. and strongly-interacting quantum gases in optical
1 ibard and co-workers for a BEC of 3·105 atoms of
0 lattices. The behavior of a Bose-Einstein conden- 87Rb [7, 8]. The resulting axially-symmetric trap
0 sate under rotation is essential for understanding
withasmallquarticanharmonicityinthexy-plane,
1 many fundamental BEC phenomena. For instance,
seen by individual atoms, has the form
:
v its response to rotation represents one of the sem-
Xi inal hallmarks of superfluidity [5]. However, once V = MRb(ω2 −Ω2)r2 +MRbω2z2+kr4 , (1)
a harmonically trapped Bose-Einstein condensate BEC 2 ⊥ ⊥ 2 z 4 ⊥
r
is rotated critically, i.e. the rotation frequency be-
a
with the perpendicular radius r = x2+y2, as
comessolargethatitcompensatestheradiallycon- ⊥
well as the trap frequencies ω = 2π × 64.8 Hz,
fining harmonic frequency, the system turns out to ⊥ p
ω = 2π × 11.0 Hz, and the trap anharmonic-
beradiallynolongertrapped. Intheabsenceofad- z
ity k = k = 2.6 × 10−11 Jm−4. Further-
ditionalpotentialtermsthecondensatewouldstart BEC
more, the rotation frequency Ω, which is measured
in units of ω , i.e. it is expressed by the ratio
∗Correspondingauthor. ⊥
r = Ω/ω , represents the tunable control param-
Email address: [email protected] (AntunBalaˇz) ⊥
1URL:http://www.scl.rs/ eter which could be experimentally increased up
Preprint submitted toPhysics LettersA January 11, 2010
to r = 1.04. This highest possible rotation fre- strate that this path-integral approach allows to
quency seems to coincide with an instability which determine a huge number of both eigenvalues and
followsfromaThomas-FermisolutionoftheGross- eigenfunctions of quantum systems with very high
Pitaevskii equation [9]. accuracy using exact numerical diagonalization.
As long as we can ignore the presence of two- In this Letter we show how this path-integral
particle interactions and approximately describe approach is applied for studying both global and
the system with the ideal Bose gas, its many- local properties of fast-rotating Bose-Einstein con-
particlepropertiesinthe grand-canonicalensemble densates. To this end we proceed as follows: Sec. 2
are exclusively derivable from one-particle states. briefly reviews the main ingredients of the path-
When considering the thermodynamic limit, usu- integral approach in the general context of ideal
ally the semiclassical approximation is applied, Bose-Einstein condensates. Then we calculate in
where the one-particle ground state E is retained Sec. 3 a large number of energy eigenvalues and
0
and treated quantum mechanically, while all one- eigenfunctions for the anharmonic one-particle po-
particlesstatesaboveE areapproximatelytreated tential(1). Afterwards,Sec.4discusseshowafinite
0
as a continuum [10]. This semiclassical approxi- number of numerically availableenergyeigenvalues
mationremains reasonablegoodirrespective of the affects the results and how they can be improved
rotation frequency Ω once the total particle num- byintroducingsystematicsemiclassicalcorrections.
ber N is large enough and the trap anharmonic- On the basis of this precise numerical one-particle
ity k small enough. The latter condition implies information, Sec. 5 studies global properties of a
that the underlying one-particle potential (1) has rotating condensate, for instance the condensation
a small curvature around its minimum. However, temperature T as a function of rotation frequency
c
in this context the question arises for which sys- Ωandtheground-stateoccupancyN /N asafunc-
0
temparameterssuchasemiclassicalapproximation tionoftemperature T. Finally, Sec.6 is devotedto
is not sufficient for a precise description of BEC thecalculationoflocalpropertiesofthecondensate,
phenomena,as wellas whenit finally breaksdown, such as density profiles and time-of-flight absorp-
requiring a full quantum mechanical treatment of tion pictures. Sec. 7 briefly summarizes the main
the system. results presented in this Letter.
In order to analyze this fundamental problem
more quantitatively, it is mandatory to deter-
2. Many- Versus One-Particle Physics
minetheone-particleenergyeigenvaluesandeigen-
functions fully quantum mechanically. To this At first we demonstrate that a precise numerical
end we apply a recently developed ultra-fast con- access to one-particle eigenstates allows the calcu-
verging path-integral approach for calculating the lation of the condensation temperature and other
imaginary-timeevolutionamplitudeofgeneralnon- thermodynamic properties of an ideal Bose gas.
relativistic many-body systems [11, 12, 13, 14, 15]. Afterwards, we briefly review details of the path-
In this approach a hierarchy of discretized effec- integral effective action approach for a numerical
tiveactionsisintroducedwherehigherterms,com- study of a 87Rb BEC in the anharmonic trap (1).
pared to the naive action, substantially reduce er-
2.1. Ideal Bose Gas
rorsincalculationsoftransitionamplitudes. Inpar-
ticular, this allows a systematic derivation of dis- For high temperatures, the grand-canonical par-
cretized effective actions of level p, which lead to a tition function of an ideal Bose gas is given by
1/Np-convergence of discretized time-sliced transi-
Z = e−β(Eν−µNν), (2)
tionamplitudeswithinthecontinuumlimitN →∞
ν
of infinitely-many time slices. In addition, the im- X
where ν enumerates all possible configurations of
provedconvergenceoftransitionamplitudes,calcu-
the system, β =1/k T represents the inverse tem-
latedwiththelevelpeffectiveaction,canbeusedto B
perature,andµ denotes the chemicalpotential. As
construct higher-order analytic approximations for
the idealbosonsdo notinteract, the systemenergy
short-time propagators. Furthermore, the N = 1
E can be expressed in terms of single-particle en-
time-slice approximation turns out to be valid for ν
ergy eigenvalues
shortimaginarytimesandusefulforsmallvaluesof
theinversetemperatureinapplicationsinquantum
E = N E , (3)
ν ν(n) n
statisticalphysics. TherecentRefs.[16,17]demon-
n
2 X
where n counts single-particle energy states, while groundstatewithintheone-particlepartitionfunc-
N =0,1,2,... and E stand for the occupancy tion, whereas the second line takes into account
ν(n) n
and the energy eigenvalue of level n, respectively. a possible macroscopic occupation of the ground
Correspondingly,thenumberofparticlesinthesys- state. The resulting total number of particles N =
tem reads −∂F/∂µ follows to be
N = N . (4)
ν ν(n) ∞
Xn N =N0+ emβµ Z1(mβ)−e−mβE0 . (9)
Thus, the grand-canonical free energy F = m=1
X (cid:2) (cid:3)
−(lnZ)/β results to be
This particle number equation servesdifferent pur-
poses in the respective phases. Within the gas
1
F = ln 1−e−β(En−µ) , (5) phase, where the macroscopic occupation of the
β
Xn h i ground state vanishes, i.e. we have N0 = 0,
Eq. (9) determines the temperature dependence of
where a subsequent Taylor expansion of the loga-
thechemicalpotentialµ. Ontheotherhand,within
rithm yields
the BEC phase the chemical potential µ coincides
1 ∞ emβµ withitsminimalvalue,i.e.theground-stateenergy
F =− Z (mβ), (6) E , so Eq. (9) yields the temperature dependence
β m 1 0
mX=1 of N0. Therefore, the value of βc =1/kBTc, which
characterizes the boundary between both phases,
usually denoted as the cumulant expansion. Thus,
followsfromEq.(9)bysettingN =0andµ=E :
0 0
the many-body thermodynamic potential (6) of an
ideal Bose gas is exclusively determined by single- ∞
particle states via the one-particle partition func- N = emβcE0Z1(mβc)−1 . (10)
tion: mX=1(cid:2) (cid:3)
Z1(β)= e−βEn. (7) We conclude that, for a given number N of ideal
n bosons, the condensation temperature can be ex-
X
Inprinciple,theaboveoutlinedexactcalculation actly calculated only if both the single-particle
of the many-body free energy allows a full numeri- ground-state energy E0 and the full temperature
cal description of ideal Bose gases, and can be also dependence of the one-particle partition function
applied for studies of dilute Bose gases in the case (7) are known.
when interactions are negligible. However, it be-
2.2. Path-Integral Approach
comes numerically very involved even for simple
trapping potentials at low temperatures. In addi- The recently developed path-integral approach
tion to this, the BEC phase transition is achieved [11, 12, 13, 14, 15] allows an efficient and fast-
onlyinthethermodynamiclimitofaninfinitenum- convergingnumerical calculation of all one-particle
ber of atoms, thus making numericalstudies of the properties of quantum systems [16, 17]. Note that
condensation increasingly difficult. Usually, this this generalnumericalapproachis suitable to treat
problem is solved by fixing the chemical potential arbitrarily-shapedtrappotentialsandisalsoappli-
µ at the low temperatures of the condensate phase cable for exact studies of many-body problems. In
to the ground-state energy, i.e. by setting µ = E , thisLetterweconsiderBECexperimentsperformed
0
andto treatthe groundstate separately,by explic- by Jean Dalibard’s group [7, 8], where a trapping
itly taking into accountits macroscopicoccupation potential with a quartic anharmonicity was real-
N . Thus, for low enough temperatures the grand- ized. For this reason,we will specialize our numer-
0
canonical free energy (6) is modified to ical approach in the following to such potentials in
order to quantitatively analyze this seminal exper-
1 ∞ emβµ iment.
F = − Z (mβ)−e−mβE0
β m 1 The quantum statistical imaginary-time transi-
mX=1 (cid:2) (cid:3) tion amplitude A(a,b;t) = hb|e−tHˆ/~|ai follows in
+N (E −µ). (8)
0 0 thecontinuumlimitN →∞ofinfinitely-manytime
In order to avoid any double-counting, we have slices
subtracted in the first line the contribution of the A(a,b;t)= lim AN(a,b;t). (11)
N→∞
3
Here the time-sliced amplitude A (a,b;t) is a is correct up to order εp, i.e. its error is propor-
N
product of short-time amplitudes corresponding to tionaltoεp+1. Ifwetakeintoaccountthepre-factor
the introduced time steps, and is expressed as a 1/(2π~ε)d/2,theaboveexpressionfortheshorttime
multiple integralof the function e−SN/~, where SN amplitude A(p) gives overall errors proportional to
isusuallycalleddiscretized(Euclidean)action[18]. εp+1/2 in d = 1. Consequently, in d = 2 spatial
For the simplest one-particle theory with the Eu- dimensions, which is the case relevant here, the er-
clidean Lagrangian L = q˙2/2 + V(q), where the rors of short-time transition amplitudes would be
coordinates are rescaled so that the mass is equal of the order εp. However, if we consider the N
to unity, the naive discretized action is given by dependence of the above expression, the errors are
alwaysproportionalto1/Np+1,sincethepre-factor
N−1 δ2 1/(2π~ε)d/2 will be absorbed in the normalization
S = n +εV(q¯ ) , (12)
N 2ε n oftheamplitudeA(a,b;t). Whenusedtocalculate
nX=0(cid:20) (cid:21) a long-time transitionamplitude, the product ofN
with the abbreviations ε = t/N, δ = q −q , short-time amplitudes will subsequently lead to er-
n n+1 n
rorsproportionalto1/Np,asdemonstratedconclu-
q¯ =(q +q )/2 andboundaryconditionsq =
n n+1 n 0
sively earlier in Refs. [11, 12, 13, 14].
a, q =b. Usingthenaivediscretizedactionleads
N
Typical values of the inverse temperature β in
to the convergence of the transition amplitudes to
BEC experiments are quite small compared to the
the continuum as slow as 1/N.
typical energy scale which is defined by the har-
Starting from the one-particle Schr¨odinger equa-
monic trap frequencies. For example, in the Paris
tion in d spatial dimensions, we have derived a hi-
erarchy of discretized effective actions S(p) with a experiment [7] the dimensionless value of ~βω⊥
N ranges between 10−3 and 10−1. Therefore, one
systematicallyimprovedconvergencetothe contin-
canimmediately use the aboveformula for the am-
uum result. The general short-time amplitude is
plitude A(p) and calculate the corresponding one-
first written in the form
particle partition function by numerically integrat-
A(q ,q ;ε) ing the diagonal amplitude A(a,a,β) over the co-
n n+1
ordinate a. For small enough β the above formula
= (2π~1ε)d/2 e−2δ~2nε−~εW(q¯n,δn;ε), (13) converges rapidly, and the amplitudes can be cal-
culated exactly for all practical purposes. We refer
were W(q¯ ,δ ;ε) stands for the ideal discretized tothisapproximationastheN =1approximation,
n n
effective action, giving exact transition amplitudes sinceitcorrespondstohavingonlyonetime-slicein
in any discretization. In fact, the above equation the standard definition of the transition amplitude
is valid for any propagation time t, not just for in the path-integral formalism.
short times, and if we were able to calculate ex- However,the directuseofthis approachhassev-
actlytheidealeffectivepotentialW,wecoulduseit eraldisadvantages. Firstofall,onestillhas to per-
directlytocalculatethelong-timetransitionampli- form an integral over the diagonal coordinate a in
tude A(a,b;t). This is practicallypossible only for order to calculate the partition function. Second,
a limited number of exactly solvable models. How- this has to be done repeatedlyfor eachvalue ofthe
ever,wehavederivedarecursiveapproachwhichal- inverse temperature β. And most importantly, one
lowsforanefficientanalyticalcalculationofahigh- has also to extract the value of the ground-state
ordershort-timeexpansionoftheeffectivepotential energyin view ofthe particlenumber equation(9).
[15]. If we insert the expansion W(p−1) calculated In principle, this is done by studying the high-β
to order εp−1 into Eq. (13), we get the short-time regime, where the short-time expansion (14) is not
amplitude valid. Althoughthisprocedureworksalsoforlower
values of β [19], it requires the numerical calcu-
A(p)(qn,qn+1;ε) lation of the one-particle partition function and a
= (2π~1ε)d/2 e−2δ~2nε−~εW(p−1)(q¯n,δn;ε), (14) dpeetraaitluerdesitnudoyrdoefrittosdoebptaenindethnecegoronutnhde-sintavteerseenteermgy-
with sufficient precision. For this reason, the al-
whichisdesignatedbypduetothefactthatW(p−1) gorithmbecomes numericallycomplex anddifficult
is multiplied by an additional factor of ε in the ex- touse,especiallyincaseswherethe groundstateis
ponent. Thisyieldsaresultfortheexponentwhich degenerate.
4
The approach based on a numerical diagonal- En/~ω⊥, k =kBEC
ization of the space-discretized evolution operator n r =0 r =1
[16, 17] effectively resolves all of the above issues. 0 1.0009731351803 0.1162667164134
Itallowstheprecisenumericalcalculationofalarge 1 2.0029165834022 0.2674689968905
numberofenergyeigenvaluesandeigenstatesofthe 2 2.0029165834022 0.2674689968905
single-particleHamiltonian. Oncecalculatedinthe 3 3.0058275442161 0.4426927375269
low-β regime, these data can be used to obtain the
4 3.0058275442161 0.4426927375270
one-particle partition function for any value of the
5 3.0067964582067 0.4725275724941
inversetemperature β in a numericallyinexpensive
6 4.0097032385903 0.6368178804983
way. Calculations based on this approach do not
7 4.0097032385903 0.6368178804984
have restrictions in the high-β regime, where, in
8 4.0116368851078 0.6848142470356
fact, they turn out to yield results with even bet-
9 4.0116368851078 0.6848142470357
ter precision. In addition, the one-particle energy
eigenfunctions obtained by the exact diagonaliza- E /~ω , k =103 k
n ⊥ BEC
tion will allow us to calculate local properties of
n r =0 r =1
Bose-Einsteincondensateswithveryhighaccuracy.
0 1.468486725893 1.162667164134
1 3.213056378201 2.674689968905
3. Numerical Calculation of Energy Eigen- 2 3.213056378201 2.674689968905
values and Eigenstates
3 5.163819069871 4.426927375269
4 5.163819069871 4.426927375270
The most efficient method for calculating prop-
5 5.406908088225 4.725275724941
erties of few-body quantum systems is the di-
6 7.282930987460 6.368178804982
rect diagonalization of the space-discretized prop-
agator in imaginary time Uˆ(t) = exp(−tHˆ/~) 7 7.282930987460 6.368178804982
8 7.690584058915 6.848142470357
[16,17,20,21,22,23]. Heretrepresentstheappro-
9 7.690584058915 6.848142470357
priatelychosenpropagationtime, for whichwe can
calculatethematrixelementsoftheevolutionoper-
atoranalyticallywithin ourpath-integralapproach
Table 1: Lowest energy levels of the xy-part of the BEC
in the N =1 approximation. potential (1) for non-rotating (r = 0) and critically ro-
On a given real-space grid defined by a set of tating (r = 1) condensate with the quartic anharmonicity
coordinates x = j∆, where ∆ denotes the spac- k = kBEC (top) and k = 103 kBEC (bottom). They are
j obtained byusinglevelp=21effective actionwiththedis-
ing, j ∈ [−L/∆,L/∆]d is a d-dimensional vector
cretizationparametersofTable3. Thespacing∆wasalways
of integers, and L denotes the spatial cutoff, ma- chosen so that L/∆ = 100, and the propagation time was
trix elements of the propagator are just short-time t=0.2fork=kBEC andt=0.05fork=103 kBEC. Errors
aregivenbytheprecisionofthelastdigit,typically10−12to
transition amplitudes 10−13,andareestimatedbycomparingthenumericalresults
obtainedwithdifferentdiscretizationparameters.
U (t)=hi∆|Uˆ(t)|j∆i. (15)
ij
Note that throughout the paper we use dimension-
lessunits,inwhichanyenergyisexpressedinterms Anothersourceoferrorsisrelatedtothefactthat
of ~ω , while the length unit is the corresponding for non-trivial potentials it is not possible to cal-
⊥
harmonic oscillator length ~/M ω . culate exactly evolution operator matrix elements.
Rb ⊥
The eigenvectors ψ (j∆) of such a discretized Higher-order effective actions minimize the errors
n p
evolutionoperatormatrix correspondto the space- for a given propagation time t and allow an opti-
discretized eigenfunctions of the original Hamilto- mal choice of this parameter as well by taking into
nian, and the corresponding eigenvalues E are re- account other discretization errors, as was shown
n
latedtotheeigenvaluesofthesingle-particleHamil- in Ref. [17]. This approach is able to give very ac-
tonianbye−tEn/~. Adetailedanalysisofdiscretiza- curate energy eigenvalues evenfor moderate values
tion errors of this method is given in Ref. [16]. It of the propagation time t of the order 0.1. Table 1
providesa practicalalgorithmforchoosingoptimal presents the first severalenergy eigenvalues for the
discretization parameters which minimizes the as- two-dimensional(xy-)partoftheBECpotential(1)
sociated errors. for the non-rotating case (r =0), as well as for the
5
E /~ω , r =1.04
n ⊥ 350
n k 103 k SC approx.
BEC BEC 300 L = 22.3, L/∆ = 100
0 -0.6617041825660 1.135693826206
1 -0.6465857464220 2.628129903790 250
2 -0.6465857464220 2.628129903790
200
3 -0.6032113415949 4.363876633929 E)
4 -0.6032113415948 4.363876633929 ρ( 150
5 -0.5349860004310 4.667653582963
100
6 -0.5349860004309 6.290444734007
7 -0.4451224795419 6.290444734007 50
8 -0.4451224795419 6.777210773169
0
9 -0.3362724309903 6.777210773169
0 20 40 60 80 100 120 140
E / −hω⊥
Table2: Lowestenergylevelsofthexy-partoftheBECpo-
Figure1: Numericallycalculateddensityofstatesforxy-part
tential (1)for over-criticallyrotating (r =1.04) condensate
withthequarticanharmonicityk=kBECandk=103kBEC ocofntdheenBsaEtCe,poobtteanitnieadlwbyithuskin=gkleBvEeCl pfo=ra21creitfficeactlliyveroatcattiionng.
accordingtothesamenumericalprocedureasinTable1.
ThediscretizationparametersareL=22.3,L/∆=100,and
t= 0.2. The dashed line is the corresponding semiclassical
approximationforthedensityofstates.
critically-rotating condensate (r = 1). The top ta-
ble givesthe energy spectrum ofthe potentialwith
theanharmonicityk =k usedintheexperiment always necessarily limited. This is easily seen from
BEC
[7], while the bottom table shows the spectrum for Fig.1,wherewecomparethedensityofstatesfora
the much larger anharmonicity k =103 k . The critically rotating condensate with the correspond-
BEC
degeneraciesof numerically obtained eigenstates in ing semiclassical approximation for the density of
all cases correspond to the expected structure of states [17]. This comparison allows to estimate
the spectrum,whichcanbe deducedfromthe sym- themaximalreliabletwo-dimensionalenergyeigen-
metry of the problem. In addition to this, the in- value E which can be obtained numerically for
max
teresting case of critical rotation (r = 1) allows a given set of discretization parameters. For ex-
a further verification of the numerical results. To ample, from Fig. 1 we can estimate E ≈ 90 for
max
this end we recall that the energy eigenvalues of a r = 1 with the discretization parameters L=22.3,
pure quartic oscillator, to which V reduces in L/∆ = 100, and t = 0.2. Table 3 gives estimates
BEC
this case, are proportional to k1/3 due to a spatial for the maximal reliable energy eigenvalue for the
rescaling in the underlying Schr¨odinger equation. anharmonicites k = k and k = 103 k for
BEC BEC
Therefore, we expect that the energy eigenvalues several values of rotation frequencies. These re-
for k =103 k are precisely 10 times largerthan sultsareobtainedfromnumericalcalculationsusing
BEC
the corresponding eigenvalues for k =k . Com- the SPEEDUP codes [24]. This table gives also an
BEC
paringthe rightmostcolumns in Table 1 we see ex- overviewoverthosediscretizationparameterswhich
actly this scaling. This demonstrates conclusively were used for a numerical diagonalization of the
that the presented method can be successfully ap- BEC potential(1)in orderto calculate bothglobal
plied also in this deeply non-perturbative parame- and local properties of the condensate throughout
ter regime. Furthermore, Table 2 gives the energy the whole paper.
spectrum of an over-critically rotating (r = 1.04) In the low-temperature limit the finiteness of
condensate,illustratingthatthesameapproachcan the number of known energy eigenstates does not
be used in this delicate regime as well. present a problem. In fact, a precise knowledge
With these results single-particle partition func- of a large number of energy eigenvalues makes this
tions Z (β) can now be calculated according to approach a preferred method for a numerically ex-
1
Eq. (7). This is especially suitable for the low- act treatment of low-temperature phenomena. On
temperatureregime,whenhigherenergylevelsgive theotherhand,thehigh-temperatureregime,where
a negligible contribution. Although the above de- thermal contributions of higher energy states play
scribed approach is able to accurately give several a significant role, is not treatable in the same way.
thousands of energy eigenvalues, their number is This regime is usually not relevant for studies of
6
k =kBEC k =103 kBEC 4. FiniteNumberofEnergyEigenvaluesand
r E /~ω L E /~ω L Semiclassical Corrections
max ⊥ max ⊥
0.0 140 14.2 190 3.90
In the previous section we have described a nu-
0.2 140 14.4 190 3.90
merical approach which is capable of providing a
0.4 140 15.0 180 3.91
large number of accurate energy eigenvalues for a
0.6 140 16.3 180 3.92
general quantum system. Its major application is
0.8 130 18.6 180 3.94
for few-body systems, where the complexity of the
1.0 90 22.3 170 3.96
algorithm is very low. For instance, we are able
1.04 90 23.2 170 3.96
tocalculatetypically104 energyeigenvaluesforthe
considered BEC potential (1). In this section we
Table 3: Maximal reliable numerically calculated energy discuss in more detail how the finiteness of nu-
eigenvalueEmax ofthexy-partoftheBECpotential(1)for merically available energy eigenstates affects the
differentvaluesofr=Ω/ω⊥,estimatedfromcomparingthe calculation of thermodynamic properties of Bose-
numerically obtained density of states ρ(E) with the semi-
Einstein condensates.
classicalapproximation. Thenumericaldiagonalization was
doneusinglevelp=21effectiveaction. Thespacing∆was As outlined in Sec. 2, the information on single-
alwayschosensothatL/∆=100,andthepropagationtime particle eigenvalues is sufficient for calculating the
was t = 0.2 for k = kBEC and t = 0.05 for k = 103 kBEC. condensation temperature according to Eq. (10).
Thetotalnumberofreliableenergyeigenstatesisinallcases
oftheorderof104. Below the condensation temperature, the ground-
state occupancy follows from solving the equation
∞
BEC experiments, but we consider it for the sake N =N + emβE0Z (mβ)−1 . (18)
0 1
of completeness. Furthermore, we want to demon-
m=1
X (cid:2) (cid:3)
stratethattheeffectiveactionapproachcanevenbe
Inpracticalcalculations,however,oneis inevitably
successfullyusedforstudyingthisparameterrange.
forced to restrict the sum over m in the cumulant
In the case that the temperature is high enough,
expansion (18) to some finite cutoff M, resulting
so that effects of higher energy eigenstates cannot
in the following approximation for the number of
be neglected, the inverse temperature β becomes
thermal atoms
a small parameter. Thus, it becomes possible to
calculate numerically the single-particle partition
M ∞
function as a sum of diagonal amplitudes, i.e. N −N = e−mβ(En−E0). (19)
0
m=1n=1
X X
Z1(β)= Ujj(β)∆d, (16) Thus the ground-state occupancy N0 depends for
j each particle number N and temperature T also
X
on this cutoff M. In particular, solving (19) for
where ∆ represents the spacing and, as before, the the (inverse) condensation temperature, obtained
values of j are defined by j ∈ [−L/∆,L/∆]d, with by demanding N0 = 0, will yield βc(M) with an
the spatial cutoff L chosen in such a way to ensure explicitdependenceonM. Theexactcondensation
the localization of the evolution matrix within the temperature βc is only obtained in the limit M →
interval [−L,L]d. Note that this can be verified, ∞.
since the evolution operator matrix elements are Fig. 2 illustrates the M-dependence resulting
calculated using the analytic approximation (14) from Eq. (19) for both a non-rotating and a crit-
with level p effective potential ically rotating condensate. As expected, the sum
saturates for high values of M to some finite num-
ber N −N . By tuning the temperature in such a
U(p)(β)= 1 e−βW(p−1)(j∆,0;~β), (17) way that th0e sum saturates at the desired value of
jj (2π~2β)d/2
the total number of atoms N in the system, which
implies N =0, one is, in principle, able to extract
0
where~β playsnow the roleofthe imaginarytime, the condensation temperature T .
c
and the influence of the kinetic term is not present Although the results in Fig. 2 suggest that this
sinceU isadiagonaltransitionamplitude(δ =0). approachcanbe applied straightforwardly,a closer
jj
7
3.0·105 3.0·105
M=30
2.9·105 2.5·105 M=1000
2.8·105 2.0·105
0 0
N - N 2.7·105 N - N 1.5·105 3.00·105
2.6·105 1.0·105 2.95·105
2.5·105 2.90·105
r=1, T=63.30 nK 0.5·105
2.4·105 r=0, T=105.18 nK 80 90 100 110 120 130 140
0
0 5 10 15 20 25 30 0 20 40 60 80 100 120 140
M Emax / −hω⊥
Figure 2: Number of thermally excited atoms N −N0 as a Figure3: NumberofthermallyexcitedatomsN−N0 calcu-
functionofthecutoffMinthecumulantexpansion(19). The latedasafunctionofthemaximalavailabletwo-dimensional
resultsaregivenforanon-rotatingcondensateatT =105.18 energy eigenvalue Emax at T = 63.30 nK. The results are
nKandforacriticallyrotatingcondensateatT =63.30nK. givenfortwodifferentvaluesofthecumulantcutoffM fora
The results are obtained by level p = 21 effective action, criticallyrotating condensate, withthe sameparameters as
andallavailablenumericaleigenstates areusedtocalculate inFig.2. Thehorizontal linecorresponds to thenumber of
N −N0. The discretization parameters were L = 14.2 for atomsN =3·105 intheexperiment[7].
r=0andL=22.3forr=1. Inbothcasesthespacingwas
chosen according to L/∆ = 100, and the propagation time
wast=0.2. A semiclassicalcorrectionto this value, canbe cal-
culated according to Ref. [10] as
lookattheresultsfornumericallycalculatedvalues ∆Z (β,E )=
1 max
ofN−N0 revealsseveralproblems that have to be d3xd3p
addressed. At first we have to investigate how the (2π~)3 e−βH(x,p)Θ(H(x,p)−Emax), (21)
resultsdependonthe numberofenergyeigenstates Z
used in the numerical calculation. Fig. 3 gives this where H(x,p) representsthe classicalHamiltonian
dependence for a critically rotating condensate at of the system, while Θ denotes the Heaviside step-
its critical temperature T =63.30 nK. We can see function.
c
thatthe dependence onthe maximalavailabletwo- For the trap potential (1) we have in z-direction
dimensional energy eigenvalue E is quite signif- a pure harmonic potential which can be treated
max
icant. The insetofthis figurerevealsanotherprob- exactly. Therefore, we focus only on the two-
lem: the value to which number N −N saturates dimensional problem in the xy-plane. In this case,
0
depends in addition on the cumulant cutoff M, as the semiclassical correction for the single-particle
explained earlier. While the M-dependence can be partitionfunction(21)canbeexpressedintermsof
dealt with by using a very large value of the cu- the complementary error function:
mulantcutoff in numericalcalculations,the depen-
denceonthemaximalenergyeigenvalueE must ∆Z(2)(β,E )= 1 e−βEmax −(1−r2)
max 1 max 2β k
beeliminatedbytakingintoaccountapropersemi- (cid:26)
classical correction to the single-particle partition (cid:2) π β(1−r2)2
+ (1−r2)2+4kEmax + e 4k
functions. kβ
Namely, the finite number of energy eigenstates p i r
β(1−r2)2
implies that the single-particle partition functions ×Erfc βEmax+ . (22)
4k
are only estimated by r !)
When this semiclassical correction is taken into
nmax
Z (β)≈ e−βEn, (20) account, the numerical results show almost no de-
1
pendenceonE ,ascanbeseenfromFig.4. Here
n=0 max
X we have used an excessively large value of the cu-
where n corresponds to the value E of the mulantcutoffM =104 inordertocompletelyelim-
max max
numerically available maximal energy eigenvalue. inate any M dependence. From the inset in this
8
4.0·105
3.0·105
With SC corr. 3.5·105
2.5·105 Without SC corr.
3.0·105
2.0·105 2.5·105
0 3.04·105 0
N - N 1.5·105 N - N 2.0·105
3.00·105 1.5·105
1.0·105 r=1.0
r=0.8
1.0·105
2.96·105 r=0.6
0.5·105
0.5·105 r=0.4
0 20 40 60 80 100 120 140 r=0.0
0 0
0 20 40 60 80 100 120 140 20 40 60 80 100 120
Emax / −hω⊥ T [nK]
Figure 4: Number of thermally excited atoms N −N0 cal- Figure 5: Number of thermally excited atoms N −N0 as
culated as a function of Emax with and without semiclas- a function of the temperature T for different values of the
sical corrections, calculated with any large cumulant cutoff rotationfrequencyandthequarticanharmonicityk=kBEC.
M = 104 to eliminate the M-dependence. The results cor- ThediscretizationparametersaregiveninTable3,andthe
respond to a critically rotating condensate with the same results are calculated by taking into account semiclassical
parameters as inFig. 2. The horizontal linecorresponds to corrections. The dashed line corresponds to the number of
N =301834 whichrepresents theexact value atTc =63.30 atoms N = 3·105 in the experiment [7]. For comparison,
nK. thefulllinesdepictthesemiclassicalresultsfromRef.[10].
graph we also see that E must be chosen in ac-
max for which the number of thermal atoms N − N
cordance with the value estimated in the previous 0
saturatespreciselyatthetotalnumberofatomsN.
section for the maximal reliable energy eigenvalue
In practice, this works the other way around: for a
obtained by numerical diagonalization. If we use a
given condensation temperature T we numerically
value E larger than this, we will be underesti- c
max calculate the particle number in the system using
mating the higher part of the energy spectra, and
Eq. (19), which gives the number of atoms in the
obtain incorrect results. For a critically rotating
system required for a condensation temperature to
condensate with the anharmonicity k = k the
BEC be equal to T . This procedure is implemented in
estimated value of E from Table 3 is around 90 c
max Fig.5forseveralvaluesoftherotationfrequencyΩ
~ω , which agrees with the results from the inset
⊥ in units of r = Ω/ω . For example, for T = 63.14
of Fig. 4. If we use this value for E and calcu- ⊥ c
max nK we see that the corresponding number of parti-
late properties of the condensate using numerically
clesisN =3·105,whichcoincideswiththevaluefor
obtainedeigenstates belowE withsemiclassical
max acriticallyrotatingcondensateintheexperimentof
correctionsaccordingtoEq.(22),wewillobtainthe
Dalibard and co-workers[7].
exact results with very high accuracy.
In principle, such a procedure is only applicable
forlow-accuracycalculationsofthecriticaltemper-
5. Global Properties of Rotating Bose- ature,sinceotherwiseonehastouseverylargeval-
Einstein Condensates ues of the cutoff M which would practically slow-
down numerical calculations. If one is interested
Inthissectionwewillapplythisapproachtocal- in more precise results, a suitable M-dependence
culate different global properties of rotating Bose- must be properly taken into account. In order to
Einstein condensates. beabletoefficientlyextractthecorrectvalueofβ ,
c
we will derive an analytical estimate of the asymp-
5.1. Condensation Temperature totic error ∆β = β −β (M) which is introduced
c c c
If we take into account semiclassical corrections by the cutoff M. Note that always ∆βc > 0, since
as explained in the previous section, we can cal- βc(M) < βc compensates the missing terms in the
culate, for instance, the condensation temperature sum (19).
ofthe condensatefor different rotationfrequencies. If we insertβ =β (M)+∆β into Eq.(19), the
c c c
This implies that we have to find the temperature error∆β canconsideredtobesmallforsufficiently
c
9
largevalueofthe cutoffM. BycomparingEq.(19)
0.05484
with the exact expression (10) we obtain
0.05482
∞ ∞ 0.05480
e−mβc(En−E0) ≈ 0.05478
m=M+1n=1 0.05476
X MX∞ M) 0.05474
∆βc m(En−E0)e−mβc(En−E0). (23) β(c 0.05472
mX=1nX=1 0.05470
The term m(E −E ) within the sum can be ob- 0.05468
n 0
tained by setting N0 = 0 and applying the partial 0.05466 kBEC, r=1.04, N=3·105
derivative ∂/∂β to Eq. (19): 0.05464
c
50 100 150 200 250 300
∂N M
−∆β ≈
c
∂β
c
Figure 6: Dependence of βc on the cumulant cutoff M for
∞ ∞
∂ an over-critically (r = 1.04) rotating condensate of N =
1−∆βc e−mβc(En−E0). (24) 3·105 atomsof87Rbwiththequarticanharmonicityofthe
∂β
(cid:20) c(cid:21)m=M+1n=1 trap k = kBEC. The discretization parameters are given
X X
in Table 3. The dashed line corresponds to a value of βc
Note that the derivative of the particle number N obtainedbyfittingthenumericalresultstothefunction(29),
with respect to β is not equal to zero, since N is whilethefulllinegivesthefittedfunctionf(M).
c
hereeffectivelydefinedbythe sum(10). Therefore,
we have instead
enoughvaluesofthe cutoffM,yieldingasasimpli-
∞ ∞ fied version of the above expression:
∂N
=− m(E −E )e−mβc(En−E0). (25)
∂βc n 0 de−(M+1)βc~ω
mX=1nX=1 ∆βc ≈−∂N/∂βc(1−e−βc~ω). (28)
Clearly, the right-hand side is a negative quantity
In order to use the derived estimates for ∆β ,
whichdoesnotdepend onM. However,itdoes de- c
apparently one would already have to know the
pendonβ andthe energyspectrumofthe system.
c
sought-after value of β as well as the difficult
Ifthesystemisclosetoad-dimensionalharmonic c
derivative ∂N/∂β . However, in practical applica-
oscillator, which is the case for the potential (1) c
tions this obstacle can be circumvented as follows.
with the small anharmonicity relevant for the ex-
Theexpressions(27)and(28)canbeusedforfitting
periment, for large values of M we have approxi-
thenumericaldataforβ (M)=β −∆β ,asisillus-
mately c c c
trated in Fig. 6. In this standard approach,all un-
∞ ∞ e−(M+1)βc~ω known values are fit parameters, obtained numer-
e−mβc(En−E0) ≈d 1−e−βc~ω , (26) ically by the least-square method. Note that not
m=XM+1nX=1 onlyβc is obtainedby suchafitting procedure,but
also other parameters, such as ∂N/∂β , or the ef-
where ω denotes an effective harmonic frequency. c
fectiveharmonicfrequencyω. Theimportantpoint
For the case of a large anharmonicity, the effec-
here is to capturethe correctM-dependence, while
tive frequency ω would depend on k, representing
all other parameters do not depend on it, so that
the harmonic expansionof the potential around its
they can be extracted by fitting. For example, in
minimum. Withsuchanestimate,Eq.(24)reduces
Fig. 6 we have used the fitting function
to
d c1e−c2(M+1)
∆βc ≈−1−e−βc~ω × f(M)=βc− 1+c3(M +1)e−c4(M+1) , (29)
e−(M+1)βc~ω which reproduces the numerical data quite accu-
. (27)
∂N/∂βc+(M +1)e−(M+1)βc~ω1−ed−~βωc~ω rately and gives high-precision results for the con-
densation temperature T . The virtue of the de-
c
The term (M + 1)e−(M+1)βc~ω in the denomina- rived estimates lies in the fact that they can be
tor of the second factor can be neglected for large usedtoextracttheinformationonthecondensation
10