Thermodynamics of a trapped Bose-condensed gas S. Giorgini1,2 , L. P. Pitaevskii2,3,4 and S. Stringari2 1 ECT∗, Villa Tambosi, Strada delle Tabarelle 286, I-38050 Villazzano, Italy arXiv:cond-mat/9704014v1 2 Apr 1997 2 Dipartimento di Fisica, Università di Trento, and Istituto Nazionale di Fisica della Materia, I-38050 Povo, Italy 3 Department 4 Kapitza of Physics, Technion, 32000 Haifa, Israel Institute for Physical Problems, 117454 Moscow, Russia Abstract We investigate the thermodynamic behaviour of a Bose gas interacting with repulsive forces and confined in a harmonic anisotropic trap. We develop the formalism of mean field theory for non uniform systems at finite temperature, based on the generalization of Bogoliubov theory for uniform gases. By employing the WKB semiclassical approximation for the excited states we derive systematic results for the temperature dependence of various thermodynamic quantities: condensate fraction, density profiles, thermal energy, specific heat and moment of inertia. Our analysis points out important differences with respect to the thermodynamic behaviour of uniform Bose gases. This is mainly the consequence of a major role played by single particle states at the boundary of the condensate. We find that the thermal depletion of the condensate is strongly enhanced by the presence of repulsive interactions and that the critical temperature is decreased with respect to the predictions of the non-interacting model. Our work points out an important scaling behaviour exhibited by the system in large N limit. Scaling permits to express all the relevant thermodynamic quantities in terms of only two parameters: the 1 reduced temperature t = T /Tc0 and the ratio between the T = 0 value of the chemical potential and the critical temperature Tc0 for Bose-Einstein condensation. Comparisons with first experimental results and ab-initio calculations are presented. 02.70.Lq, 67.40.Db Typeset using REVTEX 2 I. INTRODUCTION The experimental realization of Bose-Einstein condensation (BEC) in atomic gases [1] confined in magnetic traps has stimulated a novel interest in the theory of interacting Bose gases. In particular the recent experiments carried out at Mit on 87 23 Na [2] and at Jila on Rb [3] have provided first quantitative results for relevant thermodynamic quantities, such as the temperature dependence of the condensate fraction and the release energy. Based on the pioneering works by Bogoliubov [4] and Huang, Yang, Luttinger and Lee [5], the theory of interacting Bose gases has so far been mainly focused on uniform systems. The possibility of extending it to non-uniform gases, first explored in the works by Gross, Pitaevskii and Fetter [6,7], and stimulated by the recent experimental activity, is expected to open new challenging perspectives in many-body physics, since these systems are characterized by new length scales and symmetries. The thermodynamic behaviour of trapped Bose gases has been already the object of several theoretical works, starting from the investigation of the ideal Bose gas [8], the development of the Hartree-Fock formalism [9], the inclusion of interaction effects using the local density approximation [10,11], up to the most recent approaches based on self-consistent mean-field theory [12–15] and numerical simulations [16]. The purpose of this paper is to provide a systematic discussion of the properties of these systems at thermodynamic equilibrium. We will always assume that the dilute gas condition, na3 << 1 , (1) is satisfied. This condition is reached in all the available experiments. Here n is the atomic density whose typical values range from 1011 to 1015 cm−3 , while the triplet spin s-wave scattering length a can reach a few tens Angstrom. The corresponding value of na3 is then always extremely small (10−4 - 10−5 ) and one would naively expect that the role of interactions is a minor effect in such systems, at least for the equilibrium properties. This is not the case for these trapped atoms, because of the combined effect of trapping and BoseEinstein condensation. To illustrate this important feature let us consider the simplest case 3 of N atoms occupying the ground state of an harmonic oscillator trap. The ratio between the interaction energy and the quantum oscillator energy h̄ω behaves as: E2−body a ∼N , Nh̄ω aHO (2) where aHO = (h̄/mω)1/2 is the harmonic oscillator length. In the available experiments a/aHO is ∼ 10−3 while N ranges between 103 and 107 so that the ratio (2) can be very large, revealing that at least the ground state properties can never be calculated treating the interaction potential in a perturbative way. Differently from what happens in a homogeneous Bose gas where the density is kept fixed, the repulsion has here the effect of expanding the condensate wave function. Actually, for sufficiently large values of N, it happens that the T = 0 equilibrium properties of these new systems are governed by the competition between the 2-body force and the external potential, the kinetic energy playing a minor role. In this regime the ground state density is well approximated by the Thomas-Fermi (TF) approximation n(r) = m (µ0 − Vext (r)) , 4πh̄2 a (3) if Vext (r) ≤ µ0 and n(r) = 0 elsewhere. Here µ0 is the T = 0 chemical potential, fixed by the normalization condition. Assuming a harmonic potential Vext = m(ωx2 x2 + ωy2y 2 + ωz2 z 2 )/2, the chemical potential takes the useful expression [17,18] h̄ω a µ0 = 15N 2 aHO  2/5 , (4) where the frequency ω is the geometrical average ω = (ωx ωy ωz )1/3 (5) of the oscillator frequencies and fixes the harmonic oscillator length aHO = h̄ mω !1/2 . The root mean square (r.m.s.) radius of the condensate takes the form 4 (6) q hr 2 i  = aHO 15N a aHO 1/5 √ 2λ2 + 1 √ , 7λ2/3 (7) where λ is the anisotropy parameter of a cylindrically simmetric trap λ = ωx /ωz = ωy /ωz . In the presently available traps one has typically aHO = 10−4 cm and q hr 2 i ranges from 2 to 100 ×10−4 cm. The validity of the TF approximation (3), (4) is based on the condition Na/aHO ≫ 1 which also implies µ0 /h̄ω ≫ 1. These results have been obtained at T = 0. They play however a relevant role also at finite temperature. In particular the value of µ0 fixes an important reference in the scale of temperatures. For kB T smaller than µ0 the thermodynamic behaviour of the system is in fact affected by collective effects (phonon regime in a homogeneous superfluid). Conversely for higher temperatures it is dominated by single-particle effects. It is highly interesting to compare the value of the chemical potential µ0 with the one of the critical temperature. Due to the minor role played by interactions at high T , a reliable estimate of the critical temperature is obtained using the non-interacting model. In this model the Hamiltonian is given by H = P 2 i (pi /2m + Vext (ri )) and its thermodynamic behaviour can be worked out exactly. For sufficiently large numbers of atoms this model predicts the value [8,20] N kB Tc0 = h̄ω ζ(3) !1/3 (8) for the critical temperature. In the presently available traps one has typically h̄ω/kB = 3 −5 nK, while Tc0 ranges from 50 to 300 nK. Using the TF result (4) the ratio η between the chemical potential and the critical temperature takes the useful form η= a 2/5 µ0 = α(N 1/6 ) ≃ 2.24(a3 nT =0 (r = 0))1/6 . 0 kB Tc aHO (9) where α = 21 ζ(3)1/3 152/5 ≃ 1.57 and we have used result (3) for the density. This equation reveals an interesting behaviour of these trapped Bose gases. In fact although the gas parameter na3 is always very small, its 1/6-th power can be easily of the order of unity and in actual experiments it turns out that µ0 is 0.3-0.4 kB Tc0 . Furthermore the ratio (9) depends on the number N of atoms in the trap in an extremely smooth way (η ∼ N 1/15 ). 5 The ratio µ0 /kB Tc0 plays a crucial role to characterize the interplay between interaction and thermal effects. In particular we will show that the thermodynamic properties of the Bose condensed gases, for sufficiently large values of N, depend on the relevant parameters of the system (number N of atoms in the trap, s-wave scattering length a, oscillator length and deformation of the trap) through the dimensionless scaling parameters η. Note that such a parameter involves a different combination of the quantities N and a/aHO , as compared to the ratio (2) which is instead the natural scaling parameter for the ground state properties [17–19]. In the paper we will investigate the thermodynamic behaviour of the trapped Bose gas using a mean-field theory based on already available extensions of Bogoliubov theory to finite temperature. In particular we will work in the so called Popov approximation. The details of this approximation are discussed in Sect. 2. Here we recall some major features. At T = 0 the theory exactly accounts for the Bogoliubov spectrum of elementary excitations. At high T it approaches the finite-temperature Hartree-Fock theory [9] with which it exactly coincides above the critical temperature. These two important features (collective and single-particle effects) are the main ingredients of the Popov approximation, which is then expected to be a good approximation in the whole range of temperature except near Tc where mean-field theories are known to fail. An interesting feature of these trapped gases, which distinguishes them from uniform Bose gases, is that single-particle excitations play an important role also at low temperatures. The reason is that these systems exhibit two kinds of low lying excitations: on the one hand one has collective modes of the condensate, on the other hand one has single-particle excitations near the classical boundary. In the paper we discuss the relative importance of the corresponding contributions to thermodynamics. We find that the single-particle contribution is always relevant. In particular, even at low temperatures, the thermal energy is governed by the excitation of single-particle states differently from what happens in homogeneous Bose systems. This fact explains why in the trapped gases the Hartree-Fock approximation provides a very good description for most thermodynamic quantities in a wide range of temperatures. 6 In order to simplify the formalism and to carry out explicit calculations in a systematic way we will restrict our analysis to temperatures such that kB T >> h̄ω . (10) This assumption allows us to describe the motion of the thermally excited atoms in the WKB semi-classical approximation [21]. In a deformed trap condition (10) is separately imposed on the three oscillator energies h̄ωx , h̄ωy and h̄ωz . Current experimental investigations are carried out at temperatures well consistent with the above condition. The study of the thermodynamic behaviour for lower temperatures (of the order of the oscillator temperature) requires a quantum description of the low lying elementary modes, beyond the semi-classical approximation. These modes have been the object of recent extensive experimental and theoretical work [22–25,14]. The final equations resulting from the mean-field approach have the form of coupled equations for the order parameter and for the density of the thermally excited atoms. Their numerical solution will be systematically discussed during the paper and compared with the available experimental data as well as with the first results obtained with Path Integral Monte Carlo simulations with which we find a remarkable agreement [16]. The paper is organized as follows. Sect. 2 is devoted to the presentation of the theoretical formalism. We discuss the Popov approximation, we obtain formal results for the relevant thermodynamic quantities and we derive explicitly the scaling behaviour as well as the low T expansion of the thermodynamic functions. In Sect. 3 we discuss our results for the temperature dependence of the condensate fraction, density profile, energy, specific heat and moment of inertia. Finally in Sect. 4 we make a comparison between the properties of a trapped Bose gas and the ones of the uniform gas. II. THEORY 7 A. Self-consistent Popov approximation The self-consistent mean-field method described here follows the treatment of the inhomogeneous dilute Bose gas developed in Ref. [7] within the Bogoliubov approximation at T = 0 and its extension to the finite-temperature Hartree-Fock-Bogoliubov (HFB) approximation recently discussed by Griffin in Ref. [12]. The starting point is the grand-canonical Hamiltonian written in terms of the creation and annihilation particle field operators ψ † (r) and ψ(r) K ≡ H − µN = Z ∇2 drψ (r) − + Vext(r) − µ ψ(r) 2m ! † gZ drψ † (r)ψ † (r)ψ(r)ψ(r) . + 2 (11) In the above equation Vext (r) = m(ωx2 x2 +ωy2 y 2 +ωz2 z 2 )/2 is the anisotropic harmonic potential which provides the confinement of the atoms and g = 4πh̄2 a/m is the interaction coupling constant fixed by the s-wave scattering length a. By separating out the condensate part of the particle field operator one can write ψ(r) as the sum of the spatially varying condensate wave-function Φ(r) and the operator ψ̃(r), which acts on the non-condensate particles, ψ(r) = Φ(r) + ψ̃(r) . (12) The wave-function Φ is formally defined as the statistical average of the particle field operator, Φ(r) = hψ(r)i, and represents the order parameter of the system in the condensate phase where gauge symmetry is broken [26]. Using decomposition (12) the interaction term in (11) becomes (all quantities depend on r) ψ † ψ † ψψ = |Φ|4 + 2|Φ|2 Φ∗ ψ̃ + 2|Φ|2 Φψ̃ † + 4|Φ|2 ψ̃ † ψ̃ + Φ∗2 ψ̃ ψ̃ + Φ2 ψ̃ † ψ̃ † (13) + 2Φ∗ ψ̃ † ψ̃ ψ̃ + 2Φψ̃ † ψ̃ † ψ̃ + ψ̃ † ψ̃ † ψ̃ ψ̃ . In the following the cubic and quartic terms in ψ̃, ψ̃ † in the last line of the above equation will be treated in the mean-field approximation. For the cubic terms one gets 8 ψ̃ † (r)ψ̃(r)ψ̃(r) = 2hψ̃ † (r)ψ̃(r)iψ̃(r) ≡ 2nT (r)ψ̃(r) , ψ̃ † (r)ψ̃ † (r)ψ̃(r) = 2hψ̃ † (r)ψ̃(r)iψ̃ † (r) ≡ 2nT (r)ψ̃ † (r) , (14) where nT (r) is the non-condensate density of particles. Whereas the quartic term becomes ψ̃ † (r)ψ̃ † (r)ψ̃(r)ψ̃(r) = 4hψ̃ † (r)ψ̃(r)iψ̃ † (r)ψ̃(r) ≡ 4nT (r)ψ̃ † (r)ψ̃(r) . (15) To obtain the mean-field factorizations (14) and (15) we have neglected the terms proportional to the anomalous non-condensate density mT (r) = hψ̃(r)ψ̃(r)i, and to its complex conjugate. This approximation, which is discussed in depth in Ref. [12], corresponds to the so called Popov approximation [27], and is expected to be appropriate both at high temperatures, where nT ≫ mT , and at low temperatures where the two densities are of the same order, but both negligibly small for the very dilute systems we are considering here. Using decomposition (12) and the mean-field prescriptions (14), and (15) the grand-canonical Hamiltonian (11) can be written as the sum of four terms K = K0 + K1 + K1† + K2 . (16) The first term K0 contains only the condensate wave function Φ(r) K0 = g ∇2 + Vext (r) − µ + n0 (r) Φ(r) , drΦ (r) − 2m 2 Z ! ∗ (17) where we have introduced the condensate density n0 (r) = |Φ(r)|2 . The K1 and K1† terms are linear in the operators ψ̃, ψ̃ † K1 = Z ∇2 + Vext (r) − µ + g(n0 (r) + 2nT (r)) Φ(r) , drψ̃ (r) − 2m ! † (18) and finally K2 is quadratic in ψ̃, ψ̃ † K2 = + Z g 2 drψ̃ † (r)Lψ̃(r) + Z g 2 Z drΦ2 (r)ψ̃ † (r)ψ̃ † (r) drΦ∗2 (r)ψ̃(r)ψ̃(r) . In eq. (19), to make notations simpler, we have introduced the hermitian operator 9 (19) L=− ∇2 + Vext (r) − µ + 2gn(r) , 2m (20) n(r) = n0 (r) + nT (r) (21) where is the total particle density of the system. The K1 , K1† terms in the sum (16) vanish if Φ(r) is the solution of the equation for the condensate ∇2 + Vext (r) − µ + g(n0 (r) + 2nT (r)) Φ(r) = [L − gn0 (r)] Φ(r) = 0 , − 2m # " (22) with L defined in eq. (20), whereas the quadratic term K2 can be diagonalized by means of the Bogoliubov linear transformation for the particle operators ψ̃, ψ̃ † ψ̃(r) = X (uj (r)αj + vj∗ (r)αj† ) , j † ψ̃ (r) = X (u∗j (r)αj† + vj (r)αj ) . (23) j In the above transformation αj , and αj† are quasi-particle annihilation and creation operators which satisfy the usual Bose commutation relations, whereas the functions uj (r), vj (r) obey to the following normalization condition Z dr(u∗i (r)uj (r) − vi∗ (r)vj (r)) = δij . (24) The grand-canonical Hamiltonian (11), in the Popov mean-field approximation, takes the form KP OP OV = Z X ∇2 g drΦ(r) − ǫj αj† αj , + Vext (r) − µ + n0 (r) Φ(r) + 2m 2 j ! (25) where the real wave-function Φ(r) is the solution of eq. (22) and the quasi-particle energies ǫj are obtained from the solutions of the coupled Bogoliubov equations Luj (r) + gn0 (r)vj (r) = ǫj uj (r) Lvj (r) + gn0 (r)uj (r) = −ǫj vj (r) 10 (26) for the functions uj (r) and vj (r). In eq. (25) we have neglected the term − j ǫj P R dr|vj (r)|2 which at T = 0 accounts for the energy of the non-condensate particles in the ground state. This is a good approximation since, for dilute systems satisfying condition (1), the quantum depletion P R j dr|vj (r)|2 of the condensate is a minor effect. Explicit calculations show that this depletion is always smaller than one percent of the total number of atoms in the condensate [14]. We have therefore at T = 0, ψ(r) = Φ(r) and nT (r) = mT (r) = 0, whereas at finite temperatures we take nT and mT proportional to the thermal density of quasi-particles nT = X j (|uj |2 + |vj |2 )hαj† αj i , (27) and mT = 2 X j uj vj∗ hαj† αj i . (28) Equations (26) have the typical form of the equations of the random phase approximation and have been recently solved in Ref. [25] at T = 0; the resulting energies ǫj give the excitation frequecies of the collective modes of the system. Very recently, first solutions have become available also at finite temperature [14]. The coupled equations (26) can be solved very easily if one uses the semiclassical [13] or local density approximation [10,11]. We write uj (r) = u(p, r)eiϕ(r) vj (r) = v(p, r)eiϕ(r) , (29) where the impulse of the elementary excitation is fixed by the gradient of the phase p = h̄∇ϕ and satisfies the condition p ≫ h̄/aHO . We also assume that u(p, r), v(p, r) (and hence nT ) are smooth functions of r on length scales of the order of the harmonic oscillator length aHO = (h̄/mω)1/2 . Since the typical size of the thermally excited cloud is of the order of the classical thermal radius (2kB T /mω 2 )1/2 , the semiclassical approximation makes sense only if kB T ≫ h̄ω [21]. The sum over the quantum states j labelling the solutions of the Bogoliubov equations (26) is replaced by the integral dp/(2πh̄)3 . By neglecting derivatives 11 R of u and v and second derivatives of ϕ it is possible to rewrite the Bogoliubov-type equations (26) in the semiclassical form p2 + Vext (r) − µ + 2gn(r) u(p, r) + gn0 (r)v(p, r) = ǫ(p, r)u(p, r) 2m ! (30) 2 ! p + Vext (r) − µ + 2gn(r) v(p, r) + gn0 (r)u(p, r) = −ǫ(p, r)v(p, r) . 2m The solution of equations (30) is easily obtained and one finds the following result for the spectrum of the elementary excitations  p2 ǫ(p, r) =  + Vext (r) − µ + 2gn(r) 2m !2 1/2 − g 2 n20 (r) , (31) with the coefficients u(p, r) and v(p, r) given by p2 /2m + Vext (r) − µ + 2gn(r) + ǫ(p, r) 2ǫ(p, r) 2 p /2m + Vext (r) − µ + 2gn(r) − ǫ(p, r) v 2 (p, r) = 2ǫ(p, r) gn0(r) . u(p, r)v(p, r) = − 2ǫ(p, r) u2 (p, r) = (32) The spectrum (31) for the elementary excitations is significant only for energies ǫ(p, r) ≫ h̄ω. In fact the low-lying modes of the system are not properly described in the semiclassical approximation and their study requires the solution of the coupled Bogoliubov equations (26) [25,14]. As already mentioned above, the condition ǫ(p, r) ≫ h̄ω restricts the study of thermodynamics to temperatures much larger than the oscillator temperature h̄ω/kB . In practice the resulting analysis is useful only for systems with large values of N, where the critical temperature is much larger than the oscillator frequency (see eq.(8)). It is interesting to analyze in more detail the excitation energies (31). First of all in the absence of interparticle interactions (g = 0) the elementary excitations reduce to the usual single particle spectrum in the presence of the external potential Vext ǫ0 (p, r) = p2 + Vext (r) − µ , 2m 12 (33) where the chemical potential is given by µ = 3h̄ω̄/2, and ω̄ = (ωx + ωy + ωz )/3 is the arithmetic average of the oscillator frequencies. In the limit Na/aHO ≫ 1 (Thomas-Fermi limit) we can distinguish two regions in space depending on whether the coordinate r lies inside or outside the condensate cloud. Inside the condensate region one can neglect the kinetic term in eq. (22) and one has approximately g(n0 + 2nT ) = µ − Vext . In this region eq. (31) becomes ǫinside (p, r) = 1/2 p  2 p + 4mgn0 (r) , 2m (34) which corresponds to the local-density form of the Bogoliubov spectrum. For small p eq. (34) gives the phonon-like dispersion ǫinside (p, r) = c(r)p where c(r) = q gn0 (r)/m is the local velocity of sound. Far from the condensate cloud one instead has Vext >> 2g(n0 + nT ) and we are left with ǫoutside (p, r) = p2 + Vext (r) − µ , 2m (35) which has the same form as the independent particle spectrum (33), but now the chemical potential contains the effects of interaction and is in general much larger than its noninteracting value. Finally it is worth noting that when the energy ǫ(p, r) is much larger than the chemical potential µ one obtains the Hartree-Fock type spectrum [9] ǫHF (p, r) = p2 + Vext (r) − µ + 2gn(r) . 2m (36) In the semiclassical scheme the canonical transformation (23) becomes the familiar Bogoliubov transformation in momentum space relating the particle and quasi-particle creation and annihilation operators ψ̃(p, r) = u(p, r)α(p, r) + v(p, r)α†(−p, r) , ψ̃ † (p, r) = u(p, r)α†(p, r) + v(p, r)α(−p, r) . (37) The quasi-particle distribution function f (p, r) is obtained from the statistical average hα† (p, r)α(p, r)i, given by the usual Bose distribution 13 f (p, r) = hα† (p, r)α(p, r)i = 1 , exp(ǫ(p, r)/kB T ) − 1 (38) whereas the particle distribution function F (p, r) is given by the statistical average hψ̃ † (p, r)ψ̃(p, r)i of the particle field operators in the coordinate-momentum representation. By making use of the transformation (37) one gets   F (p, r) = hψ̃ † (p, r)ψ̃(p, r)i = u2 (p, r) + v 2 (p, r) f (p, r) ∂ǫ(p, r) =− ∂µ ! f (p, r) , (39) n0 ,nT with ∂ǫ(p, r) ∂µ ! n0 ,nT =− p2 /2m + Vext (r) − µ + 2gn(r) . ǫ(p, r) (40) It is worth noticing that the particle and quasi-particle distribution functions differ only in the low-energy regime. In fact at high energies the excitation spectrum takes the HF form (36), ∂ǫ/∂µ = −1 and hence F (p, r) = f (p, r). Conversely in the phonon regime one has F (p, r) = mc(r)f (p, r)/p. The non-condensate density nT can be obtained in a self consistent way by integrating over momenta the particle distribution (39) nT (r) = Z dp F (p, r) , (2πh̄)3 (41) and finally the chemical potential µ is fixed by the normalization condition N = N0 (T ) + NT = Z drn0 (r) + Z drnT (r) , (42) where N0 (T ) and NT are respectively the total numbers of atoms in the condensate and out of the condensate in equilibrium at temperature T . Equations (22), (31) and (38)-(42) must be solved self-consistently. The solution of these equations yields directly, for a given temperature T , the chemical potential µ, the condensate density n0 (r) and the density of thermally excited atoms nT (r). Starting from the results discussed above it is possible to study the thermodynamic properties of the system whose excited states can be classified in terms of the states of a gas 14 of undamped quasi-particles with energy spectrum (31). The entropy of the system is then readily obtained through the combinatorial formula [28] S = kB Z ǫ(p, r) f (p, r) − log(1 − exp(−ǫ(p, r)/kB T )) kB T drdp (2πh̄)3 ! . (43) One must notice that the dependence of the entropy on T is determined also through the temperature dependences of the quantities µ, n0 and nT which enter the excitation spectrum ǫ(p, r) (see eq. (31)). Starting from the entropy one can work out all the thermodynamic quantities. For example the total energy of a N-particle system at temperature T can be calculated using the thermodynamic relation E(T ) = E(T = 0) + Z 0 T dT ′ CN (T ′ ) , (44) where CN (T ) = T ∂S ∂T ! (45) N is the specific heat at constant N and E(T = 0) is the energy at zero temperature given by the well known Gross-Pitaevskii energy functional [17,19]. E(T = 0) = Z g ∇2 + Vext (r) + n0 (r) Φ(r) , drΦ(r) − 2m 2 ! (46) where Φ(r) is here the condensate wave function at T = 0. In the last part of this section we discuss the rotational properties of a trapped Bose gas which point out the effects of superfluidity [29]. The moment of inertia Θ, relative to the z direction, is defined as the linear response of the system to the perturbation field −ΩLz according to the relation hLz i = ΩΘ , (47) where Lz is the z-component of the angular momentum. In terms of the particle field operator one has Lz = Z drψ † (r) ℓz ψ(r) , 15 (48) where ℓz = −ih̄(x∂/∂y − y∂/∂x). The statistical average in eq. (47) must be calculated in the ensemble relative to the perturbed grand-canonical Hamiltonian K ′ (Ω) = K − ΩLz . In the frame rotating with angular velocity Ω the quasi-particles have energies  p2 + Vext (r) − µ + 2gn(r) ǫΩ (p, r) =  2m !2 and are distributed according to the Bose function fΩ (p, r) = − 1/2 g 2 n20 (r) − Ω(r × p)z , 1 . exp(ǫΩ (p, r)) − 1 (49) (50) In the limit of an axially symmetric trap (ωx → ωy ), the condensate wave function does not contribute to the average hLz i which is determined only by the angular momentum carried by the elementary excitations: hLz i = drdp (r × p)z fΩ (p, r) . (2πh̄)3 Z (51) To lowest order in Ω one obtains the following relevant result for the moment of inertia Θ Θ=− Z ∂f (p, r) drdp (r × p)2z . 3 (2πh̄) ∂ǫ(p, r) (52) Eq. (52) can be also written in the useful form Θ= Z dr(x2 + y 2)ρn (r) , (53) where ρn (r) = − Z dp p2 ∂f (p, r) (2πh̄)3 3 ∂ǫ(p, r) (54) provides a generalization of the most famous Landau formula for the density of the normal component to non homogeneous systems. The deviation of the moment of inertia from the rigid value Θrig = m Z dr(x2 + y 2 )n(r) , provides a ”measure” of the superfluid behaviour of the system (see section 3-C). 16 (55) B. Hartree-Fock approximation In the previous section we have seen that outside the condensate and in general for high excitation energies the Popov quasi-particle spectrum (31) reduces to the familiar HartreeFock spectrum (36). The self-consistent HF method has been already employed to study the thermodynamic properties of inhomogeneous dilute Bose systems [9,15]. The HF approximation can be easily obtained by neglecting the small ”hole” components vj in the first of the Bogoliubov equations (26). The resulting equation then takes the form ∇2 + Vext (r) − µ + 2gn(r) uj (r) = ǫj uj (r) . Luj (r) ≡ − 2m ! (56) Notice that the Hamiltonian (56) differs from the effective Hamiltonian (22) yielding the condensate wave function (for a more detailed derivation and discussion of the Hartree-Fock approximation for Bose gases see Ref. [9]). By applying the semiclassical prodedure to solve eq. (56), one finds the energy spectrum ǫHF (p, r) of eq. (36). In the HF approximation there is no difference between particles and quasi-particles and therefore the distribution functions f (p, r) and F (p, r) coincide (∂ǫHF /∂µ = −1) F (p, r) = f (p, r) = 1 . exp(ǫHF (p, r)/kB T ) − 1 (57) Results (43)-(45) for the thermodynamic quantities obtained in section 2-A hold also in the HF approximation if one simply replaces the excitation spectrum ǫ(p, r) with the corresponding HF spectrum (36). By calculating the statistical average of the grand-canonical Hamiltonian (11) over the eigenstates eq. (56) one obtains the following explicit expression for the energy of the system in the HF approximation EHF ∇2 drdp p2 = drΦ(r) − Φ(r) + F (p, r) 2m (2πh̄)3 2m Z gZ dr(2n2 (r) − n20 (r)) . + drVext(r)n(r) + 2 Z ! Z (58) In the above equation the first two terms give respectively the kinetic contribution from the condensate and non-condensate particles. The third term gives the confining energy 17 averaged over the total density n, and the last term represents the interaction contribution. The HF energy functional (58) is consistent with the thermodynamic relations ((43)-(45)). Results (52)-(54) for the moment of inertia also hold in the HF approximation with the distribution f (p, r) given by (57). In the HF approximation one easily finds ρn (r) = mnT (r) and hence the moment of inertia is given by ΘHF = m Z dr (x2 + y 2 ) nT (r) (59) consistently with the finding of Ref. [29]. The identification of the normal component ρn with the non-condensate density is a peculiar property of the HF approximation. In the Popov approach these two quantities behave differently at low temperatures because of the presence of phonon-type excitations. C. The non-interacting model The results for the thermodynamic functions of the trapped ideal gas are easily obtained by setting the coupling constant g equal to zero in the equations derived above. The equation for the condensate (22) takes the form ∇2 − + Vext (r) − µ Φ(r) = 0 , 2m ! (60) whose solution is the gaussian function Φ(r) = q N0  mω πh̄ 3/4 m exp − (ωx x2 + ωy y 2 + ωz z 2 ) 2h̄   , (61) normalized to the number of condensate atoms dr|Φ|2 = N0 , whereas the chemical potential R is µ = h̄(ωx + ωy + ωz )/2. The semiclassical spectrum (31) takes the form (33) and the noncondensate density is given by nT (r) = Z dp (exp((p2 /2m + Vext (r) − µ)/kB T ) − 1)−1 3 (2πh̄) = λ−3 T g3/2 (exp((Vext (r) − µ)/kB T )) , 18 (62) q where λT = h̄ 2π/mkB T is the thermal wavelength and gα (z) = P∞ n=1 z n /nα is the usual Bose function. The total number of non-condensate atoms is obtained by integrating nT over space NT = Z dr nT (r) = kB T h̄ω !3 g3 (exp(µ/kB T )) . (63) In the region of temperatures where the semiclassical approximation for the excited states is valid, kB T ≫ h̄ω, one can expand the g3 (z) function around z = 1 and obtain to lowest order in N −1/3 NT 3ζ(2) ω̄ −1/3 2 = t3 + N t , N 2(ζ(3))2/3 ω (64) where we have introduced the reduced temperature t = T /Tc0 , and Tc0 = (N/ζ(3))1/3h̄ω/kB is the critical temperature in the large N limit. Finite size effects, which are accounted for by the second term in eq. (64), are proportional to the ratio ω̄/ω between the arithmetic ω̄ = (ωx + ωy + ωz )/3 and the geometric averages of the oscillator frequencies and thus depend on the trap anisotropy parameter λ = ωz /ωx = ωz /ωy , having a minimum for an isotropic trap (λ = 1). These finite size effects shift the value of the critical temperature from Tc0 towards lower temperatures. By setting the right hand side of eq. (64) equal to unity one finds the following result for the shift in the critical temperature [30] ζ(2) ω̄ −1/3 ω̄ δTc0 =− N ≃ −0.73 N −1/3 . 0 2/3 Tc 2(ζ(3)) ω ω (65) The temperature dependence of the energy per particle is easily obtained in the semiclassical approximation and one has (for t < 1) E 3ζ(4) 4 3(ζ(3))1/3 ω̄ −1/3 1/3 ω̄ −1/3 3 = t + 3(ζ(3)) N t + N . NkB Tc0 ζ(3) ω 2 ω (66) It is worth noticing that the shift (65) in the critical temperature as well as the N −1/3 terms in eqs. (64), (66) exactly coincide with the results obtained to lowest order in N −1/3 from a full quantum treatment of the non-interacting model in the large N limit [30]. For N = 1000 the predictions (64), (66) are already indistinguishable from the exact result obtained by 19 summing explicitly over the excited states of the harmonic oscillator Hamiltonian. This is well illustrated in Fig. 1 where we plot the prediction N0 /N (see eq. (64)) for the temperature dependence of the condensate fraction compared to the exact quantum calculation and the large N behaviour N0 /N = 1 − t3 obtained by setting µ = 0 in eq. (63). The accuracy of the semiclassical approximation reveals that the main origin of finite size effects in the thermodynamic properties of these systems arises from the quantum effects in the equation for the condensate, which are responsible for the proper value of the chemical potential, rather than from the quantum effects in the equation for the excited states. It is also important to notice that in the large N limit all the thermodynamic properties in the non-interacting model depend only on the reduced temperature t = T /Tc0 (see eqs. (64), (66) with N → ∞). In the next section it will be shown that in the presence of interaction the thermodynamic properties of the system depend also on the additional parameter η defined in eq. (9). Finally, the ratio Θ/Θrig of the moment of inertia to its rigid value has been calculated for the non-interacting model in Ref. [29]. In terms of the reduced temperature t this result reads 1 − t3 Θ =1− , Θrig 1 − t3 + γt4 (λN)1/3 (67) where γ = 2ζ(4)/(ζ(3))4/3 ≃ 1.69. Starting from (67) one finds that in the non-interacting model the effects of finite size are much more important for Θ than for the other thermodynamic quantities. D. Thomas-Fermi limit and scaling When the number of atoms in the condensate is enough large the Thomas-Fermi approximation, which is obtained neglecting the quantum kinetic energy term in the equation for the condensate, turns out to be very accurate [17–19]. In this limit the Popov mean-field equations show an important scaling behaviour that we have recently investi20 gated in Ref. [31]. The scaling behaviour is best illustrated by introducing the dimensionless quantities p̃i = q 1/2mkB Tc0 pi , r̃i = q m/2kB Tc0 ωi ri , µ̃ = µ/kB Tc0 , ǫ̃ = ǫ/kB Tc0 , ñT = nT (2kB Tc0 /mω 2 )3/2 /N, ñ0 = n0 (2kB Tc0 /mω 2 )3/2 /N and the reduced temperature t = T /Tc0 . By neglecting the quantum kinetic energy term in eq. (22) the mean-field equations can be rewritten in the following reduced form ñ0 (r̃) = ǫ̃(p̃, r̃) =  1 µ̃ − r̃ 2 − 2g̃ñT θ(µ̃ − r̃ 2 − 2g̃ñT ) , g̃ (68) q (69) (p̃2 + r̃ 2 − µ̃ + 2g̃(ñ0 + ñT ))2 − g̃ 2 ñ20 , and ! 1 Z ∂ǫ̃ ñT (r̃) = 3 dp̃ − f (ǫ̃/t) , π ζ(3) ∂ µ̃ (70) where the renormalized coupling constant g̃ is given by  g̃ = 8πη 5/2 /15 ≃ 5.17 N 1/6 a aHO  , (71) and is hence fixed by the interaction parameter η defined in (9). In eq. (71) θ(x) is the step function and f (ǫ̃/t) = (exp(ǫ̃/t)−1)−1 is the Bose function in terms of the reduced variables. The normalization condition for the reduced densities reads: R dr̃(ñ0 + ñT ) = 1. Equations (68)-(70) exhibit the anticipated scaling behaviour in the variables t and η [32]. Starting from their solutions one can calculate the condensate and the thermal densities as well as the excitation spectrum. This gives access to all the relevant thermodynamic quantities of the system. For example the condensate fraction is given by N0 1 =1− 3 N π ζ(3) Z ! ∂ǫ̃ f (ǫ̃/t) . dr̃dp̃ − ∂ µ̃ (72) In particular the vanishing of the right hand side of eq. (72) fixes the value of the critical temperature tc = Tc /Tc0. To the lowest order in the coupling constant g̃ one finds [13] tc ≃ 1 − 0.43η 5/2 = 1 − 1.3N 1/6 a/aHO . 21 (73) The energy of the system can be calculated starting from the standard thermodynamic relation in terms of the entropy (see eq. (44)). In units of NkB Tc0 one finds: E 5 = η+ NkB Tc0 7 Z 0 t dt′ t′ ∂s , ∂t′ (74) where the entropy per particle s is given by (see eq. (43)) s(t, η) = kB 1 π 3 ζ(3) Z dr̃dp̃  ǫ̃ f (ǫ̃/t) − log(1 − exp(−ǫ̃/t)) t  . (75) Finally, the ratio Θ/Θrig of the moment of inertia to its rigid value (see eqs. (52), (55)) takes the form 1 2 Θ = 3 Θrig 3π ζ(3) t R dr̃dp̃(x̃2 + ỹ 2)p̃2 f (ǫ̃/t)(1 + f (ǫ̃/t) R . dr̃(x̃2 + ỹ 2)(ñ0 (r̃) + ñT (r̃)) (76) In the next sections we will show that the scaling behaviour is well achieved for the choice of parameters corresponding to the available experimental conditions. Scaling is consequently expected to provide a powerful tool for a systematic investigation of these trapped Bose gases. E. Elementary excitations and low temperature behaviour In this section we discuss the thermodynamic behaviour of a trapped Bose gas at low temperature (kB T ≪ µ). Though at present this regime is not easily available in experiments, its theoretical investigation is highly interesting and reveals unexpected features concerning the competition between collective and single-particle effects which are worth discussing. We will derive our results in the framework of the scaling regime discussed in section 2-D, where one assumes that the number of atoms in the trap is large enough to justify the use of the Thomas-Fermi approximation for the condensate and that the excited states can be described in the semiclassical approximation. This allows one to write the relevant thermodynamic functions in the form of integrals in coordinate and momentum space (see eqs. (72), (74)-(76)) which are well suited for an analytic expansion at low temperature. 22 This expansion corresponds to exploring the region where kB T ≪ µ. At the same time the thermal energy kB T must be larger than the typical quanta of energy characterizing the elementary excitations of the system, in order to justify the use of the semiclassical approximation. In the absence of two-body interactions the quantum of energy is fixed by the oscillator energy h̄ω. In the presence of interactions one has to distinguish between two regions. Inside the condensate the elementary excitations are of collective nature. For a spherical trap the dispersion law of these excitations can be expressed analytically in terms of the angular momentum quantum number ℓ and of the number of radial nodes according to the law [33] ǫ(n, ℓ) = h̄ω(2n2 + 2nℓ + 3n + ℓ)1/2 . (77) Their typical energy separation is still of the order of the oscillator energy h̄ω. Near the classical boundary the relevant excitations are instead of single-particle type and in first approximation can be described by the Hartree-Fock Hamiltonian (56) which in the ThomasFermi limit, and for spherical traps takes the simplified form LT F = − where R = m ∇2 + ω 2 |r 2 − R2 | , 2m 2 (78) q 2µ/mω 2 is the classical boundary of the condensate. Actually the Hartree- Fock Hamiltonian (78) gives correctly the energy of the elementary excitations only for r ≥ R (see eq. (35)). In fact for r < R it ignores Bogoliubov effects. The potential term of the Hamiltonian (78) tends to confine the low energy single-particle states near the boundary. The typical energy gaps characterizing the low energy spectrum of (78) are of the order h̄ω(R/aHO )2/3 ∼ h̄ω(µ/h̄ω)1/3 . In the large N limit this energy is larger than the oscillator energy, but still smaller than µ. As a consequence one can find a useful range of temperatures where to apply the low t expansion. In a homogeneous gas this region would correspond to the regime of collective excitations (phonons) which entirely determine the thermodynamics of the system. In a trapped Bose gas both collective and single-particle excitations taking place near the classical boundary can be important in the determination 23 of the low temperature behaviour. Actually, as it will be shown, for most thermodynamic quantities the low t behaviour is governed by single-particle excitations. Let us first consider the thermal energy of the system. At low temperature one can neglect the t-dependence in the interaction term characterizing the excitation spectrum (69) and the integral appearing in (74) can be calculated explicitly giving the result √ 4 E p̃ + 2g̃ñ0 p̃2 5 1 Z 3 3 √ d p̃d r̃ = η + NkB Tc0 7 π 3 ζ(3) r̃<√µ̃0 exp( p̃4 + 2g̃ñ0 p̃2 /t) − 1 Z 1 p̃2 + r̃ 2 − µ̃0 3 3 + 3 d p̃d r̃ , π ζ(3) r̃≥√µ̃0 exp((p̃2 + r̃ 2 − µ̃0 )/t) − 1 (79) where we have used the dimensionless variables of section 2-D with µ̃0 = µ(T = 0)/kB Tc0 = η √ which at T = 0 fixes the size of the condensate R̃ = µ̃0 (see eq. (68)). The first integral on the right hand side (r.h.s.) of eq. (79) gives the contribution to the energy from the region inside the condensate where the Bogoliubov effects take place, the second integral yields the contribution from the region outside the condensate where the spectrum is particle-like. By √ taking the linear phonon approximation ǫ̃ = p̃ 2g̃ñ0 for the excitation spectrum one finds that the first integral on the r.h.s. of (79) has a divergent behaviour near the boundary R̃ as a consequence of the vanishing of the local sound velocity. This rules out the typical t4 phonon-type behaviour for the temperature dependence of the energy. In fact a careful expansion of the integral (79) at low t shows that the relevant region contributing to the integral is near the boundary. The expansion yields the law E 5 √ = η + A ηt7/2 , 0 NkB Tc 7 (80) where A is a numerical factor. This behaviour points out the single-particle nature of the leading contribution to the thermal energy. In terms of the density of states one finds the result, valid for ǫ → 0, g(ǫ) ∼ ǫ3/2 which should be compared with the typical phonon law g(ǫ) ∼ ǫ2 . An opposite result is obtained if one investigates the low temperature behaviour of the thermal depletion (see eq. (72)) 24 NT 1 1 p̃2 + g̃ñ0 3 3 √ √ 4 = 3 d p̃d r̃ √ 4 2 N π ζ(3) r̃< µ̃0 p̃ + 2g̃ñ0 p̃ exp( p̃ + 2g̃ñ0 p̃2 /t) − 1 Z 1 1 . + 3 d3 p̃d3 r̃ √ 2 2 π ζ(3) r̃≥ µ̃0 exp((p̃ + r̃ − µ̃0 )/t) − 1 Z (81) In fact in this case the particle distribution function F (p, r), whose integral gives the thermal depletion, is enhanced at small p by the phonon contribution (see eq. (40)) which provides the leading effect in the integral (81) according to the law: NT π2 ηt2 . = √ N 3 2ζ(3) (82) Eq. (82) provides the generalization of the Ferrel law [34] holding for Bose superfluids to a trapped Bose gas. In this case the single particle contribution arising from the boundary region is a higher order effect and behaves as t5/2 . In the scaling limit one can also obtain a useful estimate for the quantum depletion of the condensate. At T = 0 the density of non-condensate atoms is given, in the semiclassical approximation, by the integral δn(r) = R dp/(2πh̄)3 v 2 (p, r), where v 2 (p, r) is given in eq. (32). By employing the Thomas-Fermi approximation one finds the result √ 2 2 (µ̃0 − r̃ 2 )3/2 θ(µ̃0 − r̃ 2 ) , δñ(r̃) = 2 3π ζ(3) (83) already discussed in Ref. [35]. It is worth noticing that the relative local quantum depletion is given by δn(r)/n(r) = 8/3(a3 n(r)/π)1/2 which coincides with the well known result for the quantum depletion of a homogeneous Bose gas [4]. By integrating (83) over space one gets the total number of quantum depleted atoms δN = N Z η3 ≃ 0.1η 3 . dr̃δñ(r̃) = √ 6 2ζ(3) (84) The above estimate of the quantum depletion agrees reasonably well with the numerical results of Ref. [14] obtained from the full solution of the Bogoliubov equations. Let us finally discuss the low temperature behaviour of the moment of inertia. The normal component (76) gets a contribution from both the collective and single-particle excitations. One can easily show that the leading contribution to the integral (76) at low temperature is given by single-particle excitations and behaves as 25 Θ ∼ t5/2 . Θrig (85) In this case the phonon contribution is a higher order effect and behaves as t4 . III. RESULTS In this section we present the numerical results obtained by solving the equations of the Popov approximation both below and above the critical temperature. In section 3-A we discuss the temperature dependence of the condensate fraction and of the peak density in the center of the trap. Other structural properties such as the temperature dependence of the density profiles and of the mean square radii of the condensed and non-condensed part of the cloud are discussed. In section 3-B we present the results for the temperature dependence of the energy and of the specif heat Finally in section 3-C we discuss the temperature dependence of the moment of inertia. In the condensed phase, T < Tc , the self-consistent procedure is schematically divided in the following steps: i) the equation for the condensate wave-function " ∇2 − + Vext(r) − µ + g(n0 (r) + 2nT (r)) Φ(r) = 0 2m # (86) is solved, using the method described in Ref. [19], for the condensate density n0 (r) and the chemical potential µ, while keeping fixed the number of particles in the condensate, N0 (T ), and the density nT (r) of thermally excited atoms. ii) The condensate density and the chemical potential found above are used to calculate the excitation energies  p2 + Vext (r) − µ + 2gn(r) ǫ(p, r) =  2m !2 1/2 − g 2 n20 (r) , (87) which yield a new non-condensate density nT through the equation (see eqs. (38), (39)) nT (r) = Z ∂ǫ(p, r) dp − 3 (2πh̄) ∂µ ! n0 ,nT 26 (exp(ǫ(p, r)/kB T ) − 1))−1 . (88) iii) From the thermal density nT the total number of non-condensed atoms is obtained by direct integration NT = Z drnT (r) , (89) and a new value for the number of atoms in the condensate is obtained from the normalization condition N = N0 (T ) + NT . (90) iv) With the new density nT and the new value for N0 (T ) steps (i)-(iii) are repeated until convergence is reached. We stop the iterative procedure when successive iterations give values for the chemical potential and the condensate fraction, N0 (T )/N, which differ by less than a few parts in 103 . In the non-condensed phase, T ≥ Tc , the condensate density n0 (r) vanishes and the self-consistent procedure becomes easier because we do not need to solve the equation (86) for the condensate. In the present work we have studied axially symmetric harmonic traps with the confining potential given by Vext (r) = m 2 2 (ωr ρ + ωz2 z 2 ) , 2 (91) where ωr is the radial frequency in the x − y plane. Clearly the calculation can be also applied to different confining potentials Vext . Before performing the calculation described above one must specify the trap parameters (frequency ωr and deformation λ = ωz /ωr ) as well as the parameters describing the system (mass m of the particles, scattering length a and total number of atoms N). For a given temperature T the above self-consistent procedure gives directly the chemical potential, µ, the number of atoms in the condensate, N0 (T ), the condensate and non-condensate densities n0 and nT . We have first checked the validity of the semiclassical approximation employed in the present work by comparing (see Fig. 2) our results for the condensate fraction with the 27 ones recently obtained in Ref. [14] by solving the full Bogoliubov equations (26) at finite temperature. The various parameter characterizing the trap and the gas are the same in the two calculations. Despite the relatively small value of N (=2000) the agreement is good also at low T , revealing that no significant errors in the analysis are introduced by the use of the semiclassical approximation for the excited states. A. Condensate fraction and density distributions In Ref. [13] we have already studied the effect of the interatomic potential on the temperature dependence of the condensate fraction and on the critical temperature. We found that a repulsive interaction among the atoms enhances the thermal depletion of the condensate with respect to the prediction of the non-interacting model for the same number of trapped atoms. The effect results in a shift of the critical temperature towards lower temperatures. In Ref. [13] we have also given an analytic expression for the shift δTc = Tc − Tc0 of the critical temperature from the prediction kB Tc0 = h̄ω(N/ζ(3))1/3 of the non-interacting model holding in the large N limit. The result, valid to lowest order in the scattering length and for large values of N, reads δTc ω̄ a ≃ −0.73 N −1/3 − 1.3 N 1/6 . 0 Tc ω aHO (92) The first term on the right hand side (r.h.s.) of the above equation is independent of the interaction and gives the first correction to Tc0 due to the finite number of atoms in the trap (see discussion in section 2-C and [30]). The second term on the r.h.s. of eq. (92) accounts for the effect of interaction; it can be either positive or negative, depending on the sign of a, increases with N and depends only on the geometric mean ω through the harmonic oscillator length aHO = (h̄/mω)1/2 . In the currently available experimental situations the shift due to interaction is the dominant effect as soon as the number of particles is of the order of 104 √ for the JILA trap (a/aHO = 7.35 × 10−3 , λ = 8) and of the order of 105 for the MIT trap (a/aHO = 2.55 × 10−3 , λ = 18/320). The shift of the critical temperature due to interaction 28 has a quite simple physical interpretation: the presence of repulsive interactions among the atoms acts in reducing the density of particles in the center of the trap. Thus lower temperatures are needed to reach there the critical density for Bose-Einstein condensation. The opposite happens in the case of attractive interactions. In Fig. 3 we present results for the temperature dependence of the condensate fraction obtained with our self-consistent mean-field method, corresponding to N = 5×104 Rb atoms in the JILA trap, and to N = 2.9 × 107 Na atoms in the MIT trap. In the same figure we also plot the predictions for the same configurations obtained switching off the two-body interaction. The enhancement of the thermal depletion due to interaction is significant over the whole temperature range and the corresponding shift of the critical temperature is also clearly visible, being in quantitative agreement with the analytic estimate (92). In Fig. 3 we also show the accuracy of the scaling behaviour discussed in section 2-D. Both configurations correspond to the value η = 0.45 of the scaling parameter. The solid curve is the scaling function obtained from (72) with η = 0.45, and the figure clearly shows that scaling is very well verified for these configurations and that hence extremely different experimental situations can give rise to the same thermodynamic behaviour. Only very close to Tc the N = 5 × 104 points exhibit small deviations from the scaling behaviour. In fact close to Tc the scaling law (72) is approached more slowly with increasing N. It is interesting to notice that finite size effects, responsible for the deviations from the scaling law 1 − t3 in the non-interacting model (see eq. (64)), are still visible for these configurations, whereas they are strongly quenched in the presence of the interaction. In Fig. 4 we present results for the condensate fraction N0 /N as a function of the reduced temperature t for three different values of the scaling parameter η, covering the presently available experimental conditions (see also Table I). In the figure the open diamonds with the error bars are the result of the Monte-Carlo simulation of Ref. [16] which correspond to the value η = 0.33 and which are in good agreement with our predictions. The dots are the experimental results of Ref. [3]. In the experiments the number of particles N varies with t, with the value of η ranging from 0.45 to 0.39. The experiments exhibit smaller 29 deviations from the non-interacting curve with respect to the ones predicted by theory. One should however keep in mind that the measured value of T is obtained from the velocity distribution of the thermal particles after expansion. The identification of this temperature with the one of the system before expansion ignores the interaction with the condensate which is expected to produce an acceleration of the thermal cloud. A preliminary estimate based on the calculation of the interaction energy between the thermally excited atoms and the condensate, shows that for T ≃ 0.5Tc0 the final kinetic energy of the thermal cloud could be about 10% larger than its value before expansion. In Fig. 5 we show the temperature dependence of the density of particles in the center of the trap for two configurations: N = 107 Na atoms in the MIT trap and N = 5000 Rb atoms in the JILA trap. The BEC transition is clearly marked by the discontinous change of slope in the two curves. In the limit of large numbers of particles the density in the center of the trap exhibits a scaling behaviour in terms of the parameters t and η when expressed in the dimensionless form discussed in section 2-D. In Fig. 6 the peak value of the reduced density ñ(r̃ = 0) = n(r = 0)(2kB Tc0 /mω 2 )3/2 /N is reported as a function of the reduced temperature t for three values of the scaling parameter η (see also Table II). At T = 0 one has ñ(r̃ = 0) = 8πη −3/2 /15 which decreases by increasing η. For T ≫ Tc0 the system approaches a classical ideal gas in the confining potential Vext where ñ(r̃ = 0) = (πt)−3/2 . Notice that in the non-interacting model the peak density does not exhibit scaling behaviour for T < Tc0 . For example at T = 0 one has ñ(r̃ = 0) = (2/π)3/2 (N/ζ(3))1/2 which diverges as N → ∞. In Figs. 7 a-d we show the density profiles of the total, condensate and non-condensate densities for four values of T /Tc0 obtained from the self-consistent solution of the equations of the Popov approximation for the configuration corresponding to N = 5000 Rb atoms in the JILA trap. The density profiles are plotted as a function of the radial coordinate. Similar profiles are obtained for the z-direction. At T /Tc0 = 0.3 (7-a) the thermal density nT is very small compared to the condensate density which almost coincides with the total density. By increasing T /Tc0 (7-b and 7-c) nT increases significantly and at T /Tc0 = 0.9 (7-c) 30 it becomes comparable with the condensate density in the center of the trap. From Figs. 7-b, 7-c it clearly appears that the thermal density is not monotonically decreasing, and reaches a flat maximum near the boundary of the condensate. This behaviour arises from the repulsive interaction with the atoms in the condensate. Above Tc0 (7-d) the condensate has disappeared and the thermal cloud has significantly spread out. In Figs. 7 a-d the scale on the vertical axis has been kept the same to give a visible evidence of the decrease of the density in the central region of the trap. The figure clearly shows the distinct behaviour of the two components of the gas: the condensate occupies the internal region of the trap while the thermal component forms a much broader cloud. In Fig. 8 we show the mean square radii of the condensate and of the thermal component as a function of T /Tc0 for the same configuration of Fig. 7. The short-dashed line refers to the mean square radius in the x-direction of the condensate atoms. By increasing T it decreases slowly and approaches a finite value at the transition point. The average radius of the thermal cloud is represented by the long-dashed line. It is always larger than the condensate radius and increases almost linearly with temperature. The solid line corresponds to average radius of the total cloud d3 rn(r)x2 /N which coincides with the condensate radius at T = 0 R and goes over to the non-condensate radius for T > Tc . In the limit of large N the radii of the atomic cloud exhibit a scaling behaviour when expressed in the reduced units of section 2-D (r̃i = q m/2kB Tc0 ωi ri ). In Figs. 9-10 we q q show how the scaling behaviour is approached, along the x ( hx̃2 i) and z-direction ( hz̃ 2 i) respectively, for the JILA (open circles) and the MIT (solid circles) configuartions. Fig. 9 shows that the scaling behaviour for q hx̃2 i is well verified for both configurations, whereas along the z-direction the configuration with the smallest number of particles (N = 5 × 104) exhibit a slight deviation from the scaling behaviour (see Fig. 10). This reveals that in the JILA trap finite size effects along the direction of stronger confinement have not completely disappeared for this value of N. In the same figures we also plot the result for the noninteracting model (η = 0). The dotted line is the result in the large N scaling limit, while the open and solid diamonds take into account finite size effects. The deviations from the 31 scaling law are much more evident in the non-interacting case. In Fig. 11 we plot the reduced mean radius of the cloud q hr̃i2 i as a function of the reduced temperature t for three values of the scaling parameter η (see also Table III). The increase of the radius due to interaction is clearly visible by comparing the three curves with the result holding in the non-interacting model for t < 1: q hr̃i2i = q E/6NkB Tc0 = q ζ(4)/2ζ(3)t2 (see eq. (66)). Notice that for spherical traps the mean square radius is proportional to the harmonic potential energy hr̃i2i = EHO /3NkB Tc0 . B. Energy and specific heat In this section we discuss the energetics of the trapped gases by presenting results for the chemical potential, total energy, release energy and specific heat. For all these quantities the scaling behaviour is well verified in the configurations realized in the experiments and we will therefore discuss our results only in the scaling limit (see section 2-D). In Fig. 12 we report the results for the chemical potential in units of kB Tc0 as a function of the reduced temperature t for three values of the scaling parameter η (see also Table IV). Notice that for t → 0 the plotted quantity coincides with η (see definition (9)). In the classical limit, T ≫ Tc0 , the dependence on the interaction parameter disappears and one finds the classical ideal gas prediction µ/kB Tc0 = t log(ζ(3)/t3). In Fig. 13 we present the results for the total energy E of the system (see also Table V). At high temperature the behaviour is given by the classical law E/NkB Tc0 = 3t. In the BEC phase the energy per particle is significantly larger than the prediction of the non-interacting model. This feature has been confirmed by the first experimental results on the temperature dependence of the release energy, which is the sum of the kinetic and interaction energy [3]. In Fig. 13 the theoretical predictions for the release energy per particle, in units of kB Tc0 , are shown together with the experimental results of Ref. [3] (see also Table VI). Though the experimental points below Tc lie well above the non-interacting curve giving evidence for interaction effects, the accuracy is not enough to make a quantitative comparison with 32 theory. Finally in Fig. 15 we plot the specific heat per particle as a function of the reduced temperature t. These curves coincide with the derivative with respect to t of the curves of Fig. 13. The effect of interaction is to round off the peak at the transition with respect to the non-interacting model. Direct experimental results on the specific heat are not yet available. An indirect estimate can be obtained by differentiating the experimental curve of the release energy (see Ref. [3]). One should however keep in mind that the release energy is roughly a factor two smaller than the total energy. C. Moment of inertia The superfluid properties of the trapped gases can be investigated by calculating the ratio between the moment of inertia of the system Θ and its rigid value Θrig (see section 2-A). In Fig. 16 our predictions for Θ/Θrig are shown as a function of the reduced temperature t. We find that, differently from the other quantities discussed above, this ratio does not exhibit a significant dependence on the scaling parameter η. The figure also shows the very different behaviour exhibited by the non-interacting model where, even for large values of N, there is a strong dependence on the number of trapped atoms and the anisotropy λ of the trap (see eq. (67)). IV. COMPARISON WITH THE HOMOGENEOUS BOSE GAS The thermodynamic behaviour of the trapped Bose gas discussed in the previous sections exhibits important differences with respect to the case of the uniform gas. The purpose of this section is to discuss these differences and to explain their physical origin. We will mainly focus on the temperature dependence of the condensate fraction where the effects are particularly visible and significant. Let us first discuss the behaviour of the ideal non-interacting gas. In the uniform case the critical temperature for the onset of Bose-Einstein condensation is given by the well 33 known expression kB Tc0 2πh̄2 = m n ζ(3/2) !2/3 (93) as a function of the density n of the gas. In the presence of harmonic confinement the critical temperature instead depends on the number of atoms and on the oscillator frequency according to the expression (8) N kB Tc0 = h̄ω ζ(3) !1/3 , (94) holding in the large N limit. By introducing the reduced temperature t = T /Tc0 one can naturally compare the thermodynamic behaviour of the two systems. In particular the temperature dependence of the condensate fraction follows the laws N0 = 1 − t3/2 , N (95) N0 = 1 − t3 , N (96) and in the uniform and in the harmonically trapped gases respectively. Different laws for the temperature dependence of the condensate are predicted to occur by changing the form of the trapping potential [8]. A major difference between the two systems concerns the relative shape of the condensate and of the thermally excited densities (see Fig. 7). In fact, differently from the uniform gas where the two components obvously overlap everywhere, in the presence of the external harmonic potential the size of the thermal component is larger than the one of the condensate. The former is in fact of the order of the classical thermal radius (2kB T /mω 2)1/2 , while the latter is fixed by the harmonic oscillator length aHO = (h̄/mω)1/2 . As a consequence, as clearly shown by Fig. 7, most of the atoms of the thermal component lie outside the condensate and form a much more dilute gas, practically for any value of T . This behaviour has important consequences when we switch on the two-body interaction. 34 The effect of the repulsive interaction is in fact twofold. One the one hand it produces an expansion of the condensate whose size can significantly increase with respect to the one of the non-interacting model (this effect is fixed by the dimensionless parameter Na/aHO as discussed in Refs. [17–19]). On the other hand at finite temperature the repulsive forces favour the depletion of the condensate. In fact atoms are energetically favoured to leave the condensate where the repulsion is strong, and to reach the non-condensate phase which is much more dilute and interaction effects are consequently less important. This is exactly the opposite of what happens in the homogenous gas where the repulsion effects are stronger in the thermal component than in the condensate. This different behaviour can be understood by considering, for example, the Hartree-Fock expression for the interaction energy (last term of eq. (58), or, more directly, the Hartree-Fock expression for the particle distribution function (we use here the Hartree-Fock formalism because the effect can be discussed in a particularly transparent way in this approximation; analogous results are obtained in the Popov approximation as proven by the numerical results reported in Fig. 17 (a)),  F (p, r) = exp((p2 /2m + Vext − µ + 2gn)/kB T ) − 1 −1 , (97) which is obtained from eqs. (57) and (36). In the homogeneous case one has Vext = 0, µ = 2gn−gn0 and the interaction effects produce a quenching of F and hence of the thermal depletion NT = dpdr/(2πh̄)3 F (p, r) with respect to the ideal gas. In the case of harmonic R confinement one sees that the relevant region of space contributing to NT lies outside the condensate where the term 2gn is negligible and interaction effects enter only through the chemical potential. In this case the repulsive interactions produce an enhancement of F with respect to the ideal gas, since the value of µ is positive for T < Tc . Notice that the value of the chemical potential is strongly affected by the interaction effects in the condensate and consequently the final effects on the thermal depletion can be sizable. In Fig. 17 (a-b) we show explicitly the temperature dependence of the condensate fraction in the uniform (a) and in the trapped (b) gases. In order to make the comparison significant, the density of the uniform gas has been fixed equal to the peak density n(r = 0) of the trapped gas at zero 35 temperature. The figure clearly shows that the effect of the interaction goes in the opposite direction in the two cases. An important consequence of the above results is that in a trapped Bose gas the repulsive interactions have the effect of lowering the critical temperature. This can be also understood by noting that the interaction produces a general expansion of the system and one has finally to deal with a more dilute gas with respect to the non-interacting gas with the same value of N. In the uniform case, where the density is kept fixed, interaction effects are known to produce instead an increase of Tc , at least for very dilute gases. The corresponding shift cannot however be calculated in mean field theory, being associated with many-body fluctuation effects typical of critical phenomena. Also the behaviour of the condensate fraction near the critical point looks different in the two cases (see Fig. 17). In the uniform gas the mean field approach gives rise to a sizable gap in the condensate fraction at Tc . This gap, which is of the order of (a3 n)1/3 , is an artifact of the theory and vanishes when fluctuation effects are taken into account near the critical point [36]. This gap is present also in the case of the trapped gas, though it has a very small effect on the condensate fraction and on other thermodynamic functions. The reason is that the size of the condensate cloud near the transition point is small and the gap in the condensate fraction turns out to be of the order of η 10 [37]. Figure 17 points out also another important difference between the uniform and trapped gases. In the latter case phonons are always important in the calculation of the condensate fraction (this is proven by the significant difference between the prediction of the HartreeFock theory and the full Popov calculation which includes phonon effects). On the other hand in the case of harmonic confinement the Hartree-Fock approximation turns out to be an excellent approximation also at relatively low temperatures. In fact for these temperatures most of the contribution to thermodynamic properties of the trapped gas arise from excitations lying close to the classical boundary where the Hartree-Fock theory is a good approximation. It is finally worth noting that the ratio between the T = 0 value of the chemical potential 36 and the critical temperature depends on the gas parameter a3 n quite differently in the two cases. In fact while in the trapped gas this ratio behaves as (a3 n(r = 0))1/6 (see eq. (9)), in the uniform gas the dependence is much stronger and is given by (a3 n)1/3 . In conclusion we have shown that a Bose gas in harmonic traps exhibits a thermodynamic behaviour quite different from the one of the uniform gas. The difference mainly originates from the different structure of the elementary excitations and from their coupling with the condensate. This should be taken into account when one compares the predictions of manybody theories obtained in uniform gases with the properties of these novel trapped systems. V. CONCLUSIONS In this paper we have presented a systematic discussion of the thermodynamic behaviour of an interacting Bose gas confined in a harmonic trap. The formalism was based on an extension of Bogoliubov theory to non uniform systems at finite temperature. By using the semiclassical WKB approximation for the excited states we have been able to investigate various thermodynamic functions for a wide range of the relevant parameters of the system. The main results emerging from our analysis are briefly summarized below: i) For large numbers of atoms in the trap the effects of two-body interactions (limited in this work to the repulsive case) can be accounted for by the scaling parameter η= a 2/5 µ0 = 1.57(N 1/6 ) , 0 kB Tc aHO (98) which gives the ratio between the T = 0 value of the chemical potential and of the critical temperature for Bose-Einstein condensation. Typical values of η in the available experiments range between 0.38 to 0.45. Physical systems with very different values of N, scattering length and oscillator frequencies are predicted to give rise to the same thermodynamic behaviour provided they correspond to the same value of η. ii) Interactions provide a significant quenching of the condensate fraction with respect to the non-interacting model and a shift of the critical temperature towards lower values. 37 This behaviour, which differs from the one exhibited by a uniform Bose gas, points out a unique and interesting feature exhibited by these trapped systems: the occurrence of elementary excitations located outside the condensate. The thermal excitation of these modes is energetically favoured by the presence of the repulsive interaction and plays a crucial role in the thermodynamics of the system. iii) Interactions are also responsible for a significant enhancement of the thermal excitation energy. Evidence for such an effect has been recently become available from first measurements of the release energy. iv) Results have been also obtained for the moment of inertia whose deviations from the rigid value point out the occurrence of superfluidity. Such deviations are significant in a rather wide range of temperatures below Tc and have been found to depend very little on the value of the interaction parameter η. v) Explicit comparison with the predictions of Path Integral Monte Carlo simulations as well as with the full solutions of the discretized Bogoliubov equations (without using the WKB approximation) have revealed that the method presented in this paper is remarkably accurate in the relevant range of temperatures. Experimental results instead are not yet enough accurate to check the theoretical predictions in a quantitative way. ACKNOLEDGMENTS We gratefully acknowledge useful discussions with Franco Dalfovo and Allan Griffin. We would also like to thank E. Zaremba for providing the numerical results of Ref. [14]. 38 REFERENCES [1] M.H. Anderson, J.R. Ensher, M.R. Matthews, C.E. Wieman and E.A. Cornell, Science 269, 198 (1995); C.C. Bradley, C.A. Sackett, J.J. Tollet and R.G. Hulet, Phys. Rev. Lett. 75, 1687 (1995); K.B. Davis, M.-O. Mewes, M.R. Andrews, N.J. van Druten, D.S. Durfee, D.M. Kurn and W. Ketterle, Phys. Rev. Lett. 75, 3969 (1995). [2] M.-O. Mewes, M.R. Andrews, N.J. van Druten, D.M. Kurn, D.S. Durfee and W. Ketterle, Phys. Rev. Lett. 77, 416 (1996). [3] J.R. Ensher, D.S. Jin, M.R. Matthews, C.E. Wieman and E.A. Cornell, Phys. Rev. Lett. 77, 4984 (1996). [4] N. Bogoliubov, J. Phys. USSR 11, 23 (1947). [5] K. Huang and C.N. Yang, Phys. Rev. 105, 767 (1957); K. Huang, C.N. Yang and J.M. Luttinger, Phys. Rev. 105, 776 (1957); T.D. Lee and C.N. Yang, Phys. Rev. 105, 1119 (1957); T.D. Lee, K. Huang and C.N. Yang, Phys. Rev. 106, 1135 (1957). [6] E.P. Gross, Nuovo Cimento 20, 454 (1961); L.P. Pitaevskii, Zh. Eksp. Teor. Fiz. 40, 646 (1961) [Sov. Phys. JETP 13, 451 (1961)]. [7] A.L. Fetter and J.D. Walecka, Quantum Theory of Many Particle Systems (McGrawHill, New York, 1971); A.L. Fetter, Ann. Phys. (N.Y.) 70, 67 (1972); A.L. Fetter, Phys. Rev. A 53, 4245 (1996). [8] S.R. de Groot, G.J. Hooyman and C.A. Ten Seldam, Proc. R. Soc. (London), A 203, 266 (1950); V. Bagnato, D.E. Pritchard and D. Kleppner, Phys. Rev. A 35, 4354 (1987). [9] V.V. Goldman, I.F. Silvera and A.J. Leggett, Phys. Rev. B 24, 2870 (1981); D.A. Huse and E.D. Siggia, J. Low Temp. Phys. 46, 137 (1982). [10] J. Oliva, Phys. Rev. B 39, 4197 (1989). [11] T.T. Chou, C.N. Yang and L.H. Yu, Phys. Rev. A 53, 4257 (1996); T.T. Chou, C.N. 39 Yang and L.H. Yu, cond-mat/9605058. [12] A. Griffin, Phys. Rev. B 53, 9341 (1996). [13] S. Giorgini, L.P. Pitaevskii and S. Stringari, Phys. Rev. A 54, 4633 (1996). [14] D.A.W. Hutchinson, E. Zaremba and A. Griffin, cond-mat/9611023. [15] A. Minguzzi, S. Conti and M.P. Tosi, Jou. of Phys. Cond. Matter, 9, L33 (1997). [16] W. Krauth, Phys. Rev. Lett. 77, 3695 (1996). [17] M. Edwards and K. Burnett, Phys. Rev. A 51, 1382 (1995); P.A. Ruprecht, M.J. Holland, K. Burnett and M. Edwards, Phys. Rev A 51, 4704 (1995). [18] G. Baym and C. Pethick, Phys. Rev. Lett. 76, 6 (1996). [19] F. Dalfovo and S. Stringari, Phys. Rev. A 53, 2477 (1996). [20] Corrections to (8) due to finite size effects have been the object of several recent papers [30]. [21] Actually the applicability of the semiclassical approximation in the region near the boundary of the condensate implies a more severe condition (see Sect. 2-E). [22] D.S. Jin, J.R. Ensher, M.R. Matthews, C.E. Wieman and E.A. Cornell, Phys. Rev. Lett. 77, 420 (1996); D.S. Jin, M.R. Matthews, J.R. Ensher, C.E. Wieman and E.A. Cornell, preprint. [23] M.-O.Mewes, M.R.Andrews, N.J. van Druten, D.M. Kurn, D.S. Durfee and W. Ketterle, Phys. Rev. Lett. 77, 988 (1996). [24] K.G. Singh and D.S.Rokhsar, Phys. Rev. Lett. 77, 1667 (1996); S. Stringari, Phys. Rev. Lett. 77, 2360 (1996). [25] M. Edwards, P.A. Ruprecht, K. Burnett, R.J. Dodd and C.W. Clark, Phys. Rev. Lett. 77, 1671 (1996). 40 [26] P.C. Hohenberg and P.C. Martin, Ann. Phys. (N.Y.) 34, 291 (1965). [27] V.N. Popov, Functional integrals and collective excitations (Cambridge University Press, Cambridge, 1987). [28] See for example L.D. Landau and E.M. Lifshitz Statistical Physics (Pergamon, Oxford, 1980) Part 1. [29] S. Stringari, Phys. Rev. Lett. 76, 1405 (1996). [30] S. Grossmann and M. Holthaus, Phys. Lett. A208, 188 (1995) and Zeit. f. Naturforsch. 50a, 323 (1995); W. Ketterle and N.J. van Druten, Phys. Rev. A 54, 656 (1996); K. Kirsten and D.J. Toms, Phys. Rev. A 54, 4188 (1996); H. Haugerud, T. Haugset and F. Ravndal, Phys. Lett. A 225, 18 (1997). [31] S. Giorgini, L.P. Pitaevskii and S. Stringari, submitted to Phys. Rev. Lett. . [32] In the limiting case η → 0 the equations can be solved analytically. For example for the condensate fraction (see eq. (72)) one finds the expansion: N0 /N = (1−t3 )−ζ(2)ηt2 (1− t3 )2/5 /ζ(3). [33] S. Stringari, Phys. Rev. Lett. 77, 2360 (1996). [34] R.A. Ferrel, N. Meynhard, H. Schmidt, F. Schwabl and P. Szepfalusy, Ann. Phys. (N.Y.) 47, 565 (1968). [35] E. Timmermans, P. Tommasini and K. Huang, cond-mat/9609234. [36] M Bijlsma and H.T.C. Stoof, Phys. Rev. A 54, 5085 (1996). [37] Notice that the same type of difficulties arise also for the space distribution ñ0 (r̃) of the condensate density in the Thomas-Fermi limit (see eq. (68). The reason is that near the boundary critical fluctuations are important not only near the transition temperature. See K. Damle, T. Senthil, S.N. Majumdar and S. Sachdev, Europhys. Lett. 36, 7 (1996). 41 FIGURES FIG. 1. Temperature dependence of the condensate fraction in the non-interacting limit. The solid circles correspond to the exact quantum calculation for N = 1000 atoms in the JILA trap and the solid line to the semiclassical approximation of eq. (64). The short-dashed line refers to the large N limit. FIG. 2. Condensate fraction as a function of T /Tc0 . The solid circles are the result of the calculation of Ref. [14] for N = 2000 Rb atoms in an isotropic trap with oscillator length aHO = 7.62×10−5 cm, obtained without using the semiclassical approximation. The solid line is the result for the same configuration employing the semiclassical approximation. Also plotted is the curve of the non-interacting model in the large N limit (short-dashed line). FIG. 3. Theoretical predictions for the condensate fraction as a function of T /Tc0 . The open √ circles correspond to N = 5 × 104 Rb atoms in the JILA trap (a/aHO = 7.35 × 10−3 , λ = 8). The solid circles correspond to N = 2.9 × 107 Na atoms in the MIT trap (a/aHO = 2.55 × 10−3 , λ = 18/320). The solid line corresponds to the scaling limit with η = 0.45. The dotted line is the 1 − t3 curve of the non-interacting model in the large N limit. The open and solid diamonds correspond to N = 5 × 104 and N = 2.9 × 107 non-interacting particles in the JILA and MIT traps respectively. FIG. 4. Condensate fraction in the scaling limit as a function of T /Tc0 for three values of η. Solid line: η = 0.45, long-dashed line: η = 0.39, short-dashed line: η = 0.31. Open diamonds: PIMC results of Ref. [16]. Solid circles: experimental results from Ref. [3]. The dotted line refers to the non-interacting model in the large N limit. FIG. 5. Peak density as a function of T /Tc0 . Solid line: N = 107 Na atoms in the MIT trap, long-dashed line: N = 5 × 103 Rb atoms in the JILA trap. FIG. 6. Reduced peak density as a function of T /Tc0 for three values of the scaling parameter η (see Fig. 4). 42 FIG. 7. Density profiles as a function of the radial coordinate at z = 0 for N = 5×103 Rb atoms in the JILA trap at different temperatures. Solid line: total density, long-dashed line: condensate density, dotted line: non-condensate density. Lengths are in units of the radial oscillator length ar = (h̄/mωr )1/2 . Densities are in units of a−3 r . FIG. 8. Root mean square (r.m.s.) radii, in units of ar , in the x-direction for N = 5 × 103 Rb atoms in the JILA trap as a function of T /Tc0 . Solid line: average radius of the total cloud; long-dashed line: average radius of the thermal cloud; short-dashed line: average radius of the condensate cloud. FIG. 9. r.m.s. radii in the x-direction in reduced units as a function of T /Tc0 for η = 0.45 (solid line). The open circles refer to N = 5×104 Rb atoms in the JILA trap. The solid circles correspond to N = 2.9 × 107 Na atoms in the MIT trap. The dotted line refers to the non-interacting model in the large N limit. The open and solid diamonds correspond to N = 5 × 104 and N = 2.9 × 107 non-interacting particles in the JILA and MIT traps respectively. FIG. 10. r.m.s. radii in the z-direction (see Fig. 9). FIG. 11. r.m.s. radii in reduced units as a function of T /Tc0 for three values of the scaling parameter η (see Fig. 4). The dotted line refers to the non-interacting model in the large N limit. FIG. 12. Chemical potential as a function of T /Tc0 for three values of the scaling parameter η (see Fig. 4). The dotted line refers to the non-interacting model in the large N limit. FIG. 13. Total energy of the system as a function of T /Tc0 for three values of the scaling parameter η (see Fig. 4). The dotted line refers to the non-interacting model in the large N limit. FIG. 14. Release energy of the system as a function of T /Tc0 for three values of the scaling parameter η (see Fig. 4). The solid circles are the experimental results of Ref. [3]. The dotted line refers to the non-interacting model in the large N limit. 43 FIG. 15. Specific heat as a function of T /Tc0 for three values of the scaling parameter η (see Fig. 4). The dotted line refers to the non-interacting model in the large N limit. FIG. 16. Ratio Θ/Θrigid as a function of T /Tc0 for three values of the scaling parameter η (see Fig. 4). The three curves coincide almost exactly and are represented by the solid line. The dotted line refers to N = 5 × 104 atoms in the JILA trap in the non-interacting model. The dot-dashed line refers to N = 2.9 × 107 atoms in the MIT trap in the non-interacting model. FIG. 17. Condensate fraction as a function of T /Tc0 . (a) refers to a homogeneous gas with the gas parameter a3 n = 6.2 × 10−5 , (b) refers to a trapped gas with the same value of the gas parameter in the center of the trap at T = 0, corresponding to the value η = 0.45. Solid line: Popov approximation; long-dashed line: Hartree-Fock approximation; dotted line: 1 − t3/2 (a) and 1 − t3 (b) curves of the non-interacting model in the homogeneous and trapped case respectively. 44 TABLES TABLE I. Condensate fraction T /Tc0 N0 /N N0 /N N0 /N η = 0.31 η = 0.39 η = 0.45 0.09 0.99 0.99 0.99 0.25 0.95 0.94 0.94 0.50 0.76 0.73 0.71 0.75 0.38 0.33 0.29 0.91 0.08 0.05 0.03 TABLE II. Peak density in reduced units T /Tc0 ñ(r̃ = 0) ñ(r̃ = 0) ñ(r̃ = 0) η = 0.31 η = 0.39 η = 0.45 0.25 3.40 2.43 1.95 0.50 3.19 2.28 1.81 0.75 2.55 1.79 1.41 1.00 0.30 0.28 0.26 1.25 0.15 0.15 0.15 45 TABLE III. Root mean square radius in reduced units T /Tc0 q q hr̃i2 i hr̃i2 i q hr̃i2 i η = 0.31 η = 0.39 η = 0.45 0.25 0.22 0.24 0.26 0.50 0.29 0.31 0.33 0.75 0.46 0.48 0.50 1.00 0.67 0.68 0.68 1.25 0.77 0.78 0.78 TABLE IV. Chemical potential µ/kB Tc0 µ/kB Tc0 µ/kB Tc0 η = 0.31 η = 0.39 η = 0.45 0.25 0.31 0.39 0.45 0.50 0.30 0.37 0.43 0.75 0.24 0.31 0.36 1.00 0.02 0.03 0.04 1.25 -0.70 -0.69 -0.69 T /Tc0 TABLE V. Total energy T /Tc0 E/N kB Tc0 E/N kB Tc0 E/N kB Tc0 η = 0.31 η = 0.39 η = 0.45 0.25 0.24 0.30 0.35 0.50 0.47 0.55 0.61 0.75 1.27 1.38 1.47 1.00 2.73 2.75 2.77 1.25 5.59 3.61 3.62 46 TABLE VI. Release energy T /Tc0 ER /N kB Tc0 ER /N kB Tc0 ER /N kB Tc0 η = 0.31 η = 0.39 η = 0.45 0.25 0.10 0.12 0.14 0.50 0.22 0.25 0.28 0.75 0.63 0.68 0.72 1.00 1.36 1.37 1.38 1.25 1.79 1.80 1.81 47