# Thermodynamics of a trapped Bose condensate with negative scattering length

## 1 INTRODUCTION

From 1995 we have experimental results interpreted as an evidence of Bose-Einstein condensation (BEC) in clouds of confined alkali-metal atoms . In these dilute vapors one can separate off the effect of BEC and study the role of the interaction. For the theoretical understanding of these experiments, an important contribution comes from Pitaevskii and Gross with their idea of a space-time dependent macroscopic wavefunction (order parameter) that describes the Bose-condensed phase . The experiments with alkali-metal atoms generally consist of three steps : a laser cooling and confinement in an external potential (a magnetic or magneto-optical trap), an evaporative cooling (the temperature reached is of the order of $`100`$ nK) and finally the analysis of the state of the system. Nowadays a dozen of experimental groups have achieved BEC by using different geometries of the confining trap and atomic species.

In this paper we perform a study about the thermodynamics of a Bose-Einstein condensate of <sup>7</sup>Li atoms, which is particularly interesting because of the attractive interatomic interaction: only a maximum number of atoms can form a condensate; beyond that the system collapses. Up to now there are only few experimental data supporting BEC for <sup>7</sup>Li vapors with a limited condensate number. Because of the small number of condensed atoms, the size of the Bose condensate is of the same order of the resolution of the optical images. So it is difficult to obtain precise quantitative information, like the number of condensed atoms and the BEC transition temperature (see Bradley et al. ). A detailed theoretical study of the thermodynamical properties of such systems may give useful information for future experiments with a better optical resolution.

We consider an isotropic harmonic trap and concentrate our attention on the density profiles of condensate and non-condensate fraction and on the energy spectrum of elementary excitations by using both the Bogoliubov and Popov approximations. Actually, the Popov method and its quasi-classical approximation are a good starting point to estimate the contribution of the correlations in the gas for any finite value of the temperature .

In Section II we describe the theoretical model used for our calculation and discuss in detail the properties of systems with attractive interatomic interaction. In Section III we illustrate the numerical procedure and the method we used to reduce numerical errors. Finally, in Section IV we analyze the numerical results.

## 2 THEORY

The Heisenberg equation of motion for the bosonic field operator $`\widehat{\psi }(𝐫,t)`$ of a non-relativistic system of confined and interacting identical atoms is given by 

$$i\mathrm{}\frac{}{t}\widehat{\psi }(𝐫,t)=\left[\frac{\mathrm{}^2}{2m}^2+V_0(𝐫)\mu \right]\widehat{\psi }(𝐫,t)+d^3𝐫^{}\widehat{\psi }^+(𝐫^{},t)V(𝐫,𝐫^{})\widehat{\psi }(𝐫^{},t)\widehat{\psi }(𝐫,t),$$
(1)   
where $`m`$ is the mass of the atom, $`\mu `$ the chemical potential, $`V_0`$ is the confining external potential and $`V`$ is the interatomic potential. We separate out the condensate part setting

$$\widehat{\psi }(𝐫,t)=\mathrm{\Phi }(𝐫)+\widehat{\varphi }(𝐫,t)$$
(2)   
where $`\mathrm{\Phi }(𝐫)=\widehat{\psi }(𝐫,t)=\widehat{\psi }(𝐫)`$ is the order parameter (macroscopic wavefunction) of the condensate, and $`\mathrm{}`$ is the mean value in the grand-canonical ensemble. We can obtain an equation for $`\mathrm{\Phi }(𝐫)`$ by inserting (2) into (1), by using the mean-field approximation for the fluctuations operators $`\widehat{\varphi }`$ and $`\widehat{\varphi }^+`$ , and by taking the mean value. To obtain the non-condensate fraction, we subtract from (1) its mean-value and use again the mean-field approximation. The fluctuation operator can be expanded in the following way

$$\widehat{\varphi }(𝐫,t)=\underset{j}{}\left[u_j(𝐫)\widehat{a}_je^{iE_jt/\mathrm{}}+v_j^{}(𝐫)\widehat{a}_j^+e^{iE_jt/\mathrm{}}\right],$$
(3)   
where $`\widehat{a}_j`$ and $`\widehat{a}_j^+`$ are bosonic operators and the complex functions $`u_j(𝐫)`$ and $`v_j(𝐫)`$ are the wavefunctions of the so-called quasi-particle excitations of energy $`E_j`$.

When the density of the atomic cloud is such that the scattering length and the range of the interatomic interaction are less than the average interatomic distance, the true interatomic potential can be approximated by a local pseudopotential $`V(𝐫,𝐫^{})=g\delta ^3(𝐫𝐫^{})`$ , where $`g=4\pi \mathrm{}^2a_s/m`$ is the scattering amplitude of the spin triplet channel ($`a_s`$ is the s-wave scattering length). In the typical conditions of confined <sup>7</sup>Li vapors, the density is about $`\rho =10^{20}`$ m<sup>-3</sup> and $`a_s=1.410^9`$ m, i.e. $`\rho |a_s|^3<10^6`$, thus the pseudopotential gives reliable results. So we have the Hartree-Fock-Bogoliubov equation

$$\widehat{L}u_j(𝐫)+gn(𝐫)u_j(𝐫)+gm(𝐫)v_j(𝐫)=E_ju_j(𝐫),$$

$$\widehat{L}v_j(𝐫)+gn(𝐫)v_j(𝐫)+gm(𝐫)u_j(𝐫)=E_jv_j(𝐫),$$
(4)   
where $`\widehat{L}=\frac{\mathrm{}^2}{2m}^2+V_0(𝐫)\mu +gn(𝐫)`$, and the $`u_j(𝐫)`$ and $`v_j(𝐫)`$ satisfy the normalization condition

$$d^3𝐫[u_j^{}(𝐫)u_k(𝐫)v_j^{}(𝐫)v_k(𝐫)]=\delta _{jk}.$$
(5)   
We write $`n(𝐫)=\widehat{\psi }^+(𝐫)\widehat{\psi }(𝐫)=n_0(𝐫)+n_{nc}(𝐫)`$ for the total local density, $`m(𝐫)=\widehat{\psi }(𝐫)\widehat{\psi }(𝐫)=\mathrm{\Phi }(𝐫)\mathrm{\Phi }^{}(𝐫)+m_{nc}(𝐫)`$ for the total anomalous average function, $`n_0(𝐫)=N_0|\mathrm{\Phi }(𝐫)|^2`$ for the density of the condensate fraction ($`N_0`$ is the number of condensed atoms), $`n_{nc}(𝐫)=\widehat{\varphi }^+(𝐫)\widehat{\varphi }(𝐫)`$ for the local density of the non-condensate fraction, $`m_{nc}(𝐫)=\widehat{\varphi }(𝐫)\widehat{\varphi }(𝐫)`$ for the anomalous average of the non-condensate fraction.

Nevertheless, these equations remain very complex and some simplification is usually performed.

One possible approximation, called Bogoliubov standard approximation, neglects the non-condensate fraction ($`m_{nc}=n_{nc}=0`$) in the equations. This approximation is particularly useful near zero temperature when the non-condensate is very small compared with the condensate fraction; so the condensate and its elementary excitations are separately calculated. Another possible, less drastic, simplification is called Popov approximation. It neglects only the anomalous average function $`m_{nc}`$ and it is able to give a more accurate description of the spectrum and non-condensate fraction (see Ref. 7 for the case of a homogeneous gas). Popov approximation is expected to be reliable in the whole range of temperature except near $`T_{_{BEC}}`$, where mean-field theories are known to fail. This set of equations must be solved with a self-consistent procedure that involves an iterative process. The equations are

$$\left[\frac{\mathrm{}^2}{2m}^2+V_0(𝐫)\mu \right]\mathrm{\Phi }(𝐫)+g[n_0(𝐫)+2n_{nc}(𝐫)]\mathrm{\Phi }(𝐫)=0,$$
(6)   
$$\widehat{L}_Pu_j(𝐫)+gn_0(𝐫)v_j(𝐫)=E_ju_j(𝐫),$$

$$\widehat{L}_Pv_j(𝐫)+gn_0(𝐫)u_j(𝐫)=E_jv_j(𝐫),$$
(7)   
where $`\widehat{L}_P=\frac{\mathrm{}^2}{2m}^2+V_0(𝐫)\mu +2gn(𝐫)`$. These equations are supplemented by the relation fixing the total number of atoms in the system

$$N=d^3𝐫[n_0(𝐫)+n_{nc}(𝐫)],$$
(8)   
where

$$n_0(𝐫)=N_0|\mathrm{\Phi }(𝐫)|^2$$

$$n_{nc}(𝐫)=\underset{j}{}\left(|u_j(𝐫)|^2+|v_j(𝐫)|^2\right)\widehat{a}_j^+\widehat{a}_j+|v_j(𝐫)|^2$$
(9)   
with $`\widehat{a}_j^+\widehat{a}_j=(e^{E_j/k_BT}1)^1`$ the Bose factor ($`k_B`$ is the Boltzmann constant) at temperature $`T`$. The Popov approximation was first carried out in  for repulsive interaction.

A possible way to avoid the solution of the partial differential equations (6) and (7) is to use the quasi-classical formula  for the Popov level of the excitation energies given by

$$E(𝐩,𝐫)=\sqrt{\left[\frac{p^2}{2m}+V_0(𝐫)\mu +2gn(𝐫)\right]^2g^2n_0^2(𝐫)}.$$
(10)   
With the help of this quasi-classical formula, the non-condensate density (thermal depletion) reads

$$n_{nc}(𝐫)=\frac{d^3𝐩}{(2\pi \mathrm{})^3}\left(\frac{E}{\mu }\right)\left(e^{E(𝐩,𝐫)/k_BT}1\right)^1.$$
(11)   
As suggested in Ref. 8, the quasi-classical Popov approximation is accurate only if $`E>>E_0`$, where $`E_0`$ is the lowest energy level of the elementary excitations.

We study the BEC in an external harmonic potential with spherical symmetry, which is given by

$$V_0(𝐫)=\frac{1}{2}m\omega ^2r^2,$$
(12)   
where $`r=(x^2+y^2+z^2)^{1/2}`$ and $`\omega `$ is the trap frequency. This symmetry makes easier the numerical procedure because one has to solve a one-dimensional problem.

For future reference, we recall some of the results in the case of non interacting particles (i.e. $`a_s=0`$) confined in an harmonic potential. In the limit of $`k_BT>>\mathrm{}\omega `$ one can neglect the discretization of the eigenvalues (this corresponds to the limit $`N\mathrm{}`$ and $`\omega 0`$ with $`N\omega ^3`$ constant), and one gets the following result :

$$T_{_{BEC}}^0=\frac{\mathrm{}\omega }{k_B}\left(\frac{N}{\zeta (3)}\right)^{1/3}=0.94\frac{\mathrm{}\omega }{k_B}N^{1/3},$$
(13)   
where $`T_{_{BEC}}^0`$ is the critical temperature, i.e. the temperature above which there is no a macroscopic population of the ground-state. The condensate fraction as a function of the temperature is given by

$$\frac{N_0}{N}=1\left(\frac{T}{T_{_{BEC}}^0}\right)^3.$$
(14)   
The finite size correction to lowest order in $`1/N`$ of Eq. (13) is 

$$\frac{\delta T_{_{BEC}}^0}{T_{_{BEC}}^0}=\frac{\zeta (2)}{2\zeta (3)^{2/3}}N^{1/3}=0.73N^{1/3}.$$
(15)   
The finite size correction of Eq. (14) is instead :

$$\frac{N_0}{N}=1\left(\frac{T}{T_{_{BEC}}^0}\right)^3\frac{3\zeta (2)}{2[\zeta (3)]^{2/3}}\left(\frac{T}{T_{_{BEC}}^0}\right)^2N^{1/3}.$$
(16)   
The correction in the condensate fraction due to the interaction has also been estimated as 

$$\frac{\delta T_{_{BEC}}}{T_{_{BEC}}^0}=1.3\frac{a_s}{a_h}N^{1/6},$$
(17)   
where $`a_h=(\mathrm{}/m\omega )^{1/2}`$ is the characteristic length of the trap. We will compare the numerical solution of Popov equations with these analytical estimates. Of course one should keep in mind that in a finite system no sharp transition can be present and, if present, it can only be an effect of approximation.

In the case of negative scattering length, a simple variational calculation with a Gaussian trial wavefunction shows that a harmonic trap supports a Bose condensate at $`T=0`$ with a number of bosons smaller than a $`N_c=0.67(a_h/|a_s|)`$ ) while numerical calculations give $`N_c=0.57(a_h/|a_s|)`$ . In a homogeneous gas such critical number is zero.

For sake of completeness, we observe that, when a realistic nonlocal (finite range) interaction is taken into account, a new branch of Bose condensate appears for <sup>7</sup>Li at higher density , but this non locality has a small effect on the low density branch unless the trap is much narrower of those of present experiments.

## 3 NUMERICAL PROCEDURE

To determine the order parameter of the condensate we use the steepest descent method . The ground-state energy per particle and its wavefunction is obtained by minimizing the energy functional

$$_p[\mathrm{\Phi }(𝐫)]=d^3𝐫\left[\frac{\mathrm{}^2}{2m}|\mathrm{\Phi }(𝐫)|^2+V_0(𝐫)|\mathrm{\Phi }(𝐫)|^2+\frac{gN}{2}|\mathrm{\Phi }(𝐫)|^4+gn_{nc}(𝐫)|\mathrm{\Phi }(𝐫)|^2\right].$$
(18)   
Moreover, the ground-state wavefunction must satisfy the normalization condition $`d^3𝐫|\mathrm{\Phi }(𝐫)|^2=1`$. After finite-difference spatial discretization, our function $`\mathrm{\Phi }(r)`$ is represented by the vector $`\stackrel{}{\mathrm{\Phi }}`$, whose elements are the values of $`\mathrm{\Phi }(r)`$ on a grid of points of the radial axis. The functional (18) becomes a real function $`(\stackrel{}{\mathrm{\Phi }})`$. So, starting from a normalized trial vector $`\stackrel{}{\mathrm{\Phi }}^{(0)}`$, we can obtain the minimum by iterating the map

$$\stackrel{}{\mathrm{\Phi }}^{(j+1)}=\stackrel{}{\mathrm{\Phi }}^{(j)}\stackrel{}{}(\stackrel{}{\mathrm{\Phi }}^{(j)}),$$
(19)   
and stopping the process when the change from the (j)-step to (j+1)-step is less than a fixed threshold. Note that the function must be normalized at each step. We use lattices up to $`500`$ points verifying that the results do not depend on the discretization parameters.

Also the Bogoliubov and Popov elementary excitations are obtained by finite-difference discretization  of Eq. (7) (to obtain the Bogoliubov levels we put $`n_{nc}=0`$ in these equations) and we have used a grid of $`500`$ points in the radial direction. Each eigenvector contains in its entries the discretized functions $`u_j`$ and $`v_j`$. In our case, each function takes up $`500`$ elements in a array of $`1000`$ so that the eigenvalue problem reduces to the diagonalization of a $`1000\times 1000`$ real matrix .

It is necessary to perform other two operations to get our goal. First, we have to check the numerical errors in the evaluation of the eigenfunctions at the boundary of the mesh, where they must go to zero with a decreasing oscillatory behavior; second, we must impose the normalization condition (5) and select eigenvalues and eigenvectors. We check that the sequence of maxima of each absolute value of the functions $`u_j`$ and $`v_j`$ is decreasing. When these functions are very close to zero, because of round off errors, this regular oscillatory and decreasing behavior is lost and we put the eigenfunctions equal to zero until the end of the grid.

For a given temperature $`T`$, we impose an energy threshold of $`30k_BT`$, and in any case in the sum (9) are considered only eigenstates with energy smaller then a cut off ($`70\%`$ of the value of the potential in the last point of the grid). In this way we avoid to use the highest eigenvalues, those which are most sensitive to cut off errors. The symmetry of the system allows us to decompose the series (9) and to consider separately the contribution of the levels with a fixed angular quantum number $`l`$.

Our codes have been tested with solvable models, like the non-interacting case, by comparing numerical results with analytical solutions. We find that, for a good determination of the eigenvalues, it is not necessary a great number of points but one needs a large spatial range. A good choice, in our case, is a grid of $`500`$ points in a range $`0<r<35`$ (in units of the characteristic length of the trap $`a_h=3.23`$ $`\mu `$m).

## 4 RESULTS

### 4.1 Gross-Pitaevskii equation

First of all we study the ground state properties by using the equation (6) with $`n_{nc}=0`$, that is the Gross-Pitaevskii (GP) equation. This equation describes accurately the wave function of the condensate at zero temperature because the quantum depletion is negligible. We use the following numerical values for our system: m(<sup>7</sup>Li)$`=6.941a.m.u.`$, $`\omega =878`$ Hz for the frequency of the isotropic harmonic trap and $`a_s=27a_0`$ for the scattering length ($`a_0`$ is the Bohr radius). We measure the radial coordinate $`r`$ in units of the characteristic oscillator length $`a_h`$ and the energies in units of energy oscillator quantum $`\mathrm{}\omega `$.

The GP equation is solved for different values of the total number $`N`$ of atoms in the trap. The results are shown in Fig. 1, where we plot some properties of the solutions vs. $`N`$. All these curves show the existence of a critical threshold $`N_c`$. If we choose a $`N>N_c`$ we find that the solution of the GP equation is a sort of delta function centered at $`r=0`$ and the energy per particle goes to $`\mathrm{}`$. The kinetic energy cannot compensate the increasing interaction energy any more: the energy functional has in these conditions only an absolute minimum that correspond to the collapsed system. Our numerical calculations show that the threshold is equal to $`N_c=1287`$, in full agreement with a previous numerical estimation ; instead, one finds $`N_c=1514`$ by using the approximate variational formula with a Gaussian trial wavefunction.

### 4.2 Far from the critical threshold $`N_c`$

Now we analyze in detail the system with $`N=1000`$, a number of atoms sufficiently far from the critical threshold $`N_c`$. We solve self-consistently the Popov equations (6) and (7) for some value of the temperature. In this way, we obtain information about the density profiles of the condensate and non condensate fractions (9) and the spectrum of elementary excitations.

It is interesting to study the temperature dependence of the elementary excitations and the modification induced by the self-consistent Popov procedure. The changes are very small in both cases, as reported in Fig.2a and Tab.1. The energy spectrum is well represented by a linear function of $`n`$, i.e. by the spectrum of a harmonic oscillator with a small shift. A peculiar elementary excitation is that with radial quantum number $`n=0`$ and angular quantum number $`l=1`$. This is the dipole mode frequency at which the center of mass of the atomic cloud oscillates. Due to the generalized Kohn theorem such dipole mode should have the same frequency of the trap, but the Popov approximation does not fully respect this theorem because the thermal cloud is static . Except for the $`n=1`$ and $`l=0`$ mode, the Popov method gives eigenvalues lower than Bogoliubov ones. This effect is a consequence of the shape of the effective potential of the Bose condensate. In the Popov approximation, the effective potential is given by $`V_0(𝐫)+2g(n_0(𝐫)+n_{nc}(𝐫))`$, thus it depends on the non condensate density profile too. This term gives, in the case of $`l`$0, an effective potential with a larger parabolic shape and then with lower energy levels. The mode with quantum numbers $`n=1`$ and $`l=0`$, instead, is affected by the shape of the potential near $`r=0`$, where the term $`2g(n_0(𝐫)+n_{nc}(𝐫))`$ produces a significant correction to the harmonic potential. So the presence of the non-condensed density profile gives a deeper well and a more bound state. In Fig. 2b we show the lowest energy levels as a function of the temperature. One sees that by increasing the temperature they approach even better the energy levels of the harmonic oscillator. This behaviour is due to the fact that by increasing the temperature the number of condensed atoms decreases and the nonlinear term of the eq. (7) is reduced.

In Fig. 3 the density profiles of the condensed and non condensed atoms are shown. We plot in Fig. 3a the profile of the condensate fraction for some value of temperature from $`0`$ to $`58`$nK. Note that the temperature does not modify the shape of the condensate density and its radial range. This is not the case of the non condensed density profile. As shown in Fig. 3b, the increasing temperature produces a sharper maximum at a larger value of $`r`$. Finally, in Fig. 4 we show the number of non condensed atoms for each value of the angular quantum number $`l`$, i.e. the population of the mode with a fixed quantum number $`l`$. For three values of the temperature, we compare the results in Bogoliubov and Popov approximation. The contribution for $`l=0`$ in Popov approximation is smaller than in Bogoliubov approximation because of the increase in the excitation energy obtained by the Popov method. We can see a change in the shape of $`N_{nc}(l)`$ with temperature, related to the modification of the non condensed density profiles. The more visible differences from Bogoliubov to Popov are due to the modes with lower quantum number $`l`$. Also this behavior is related to the form of the effective potential because the higher energy levels are not affected to the changes of the form of the potential as the lowest levels in the spectrum.

In Tab. 2 we compare the Bogoliubov and Popov results. The differences become more pronounced as temperature is increased.

We summarize the results of our calculations with Popov method in Fig. 5, where we plot the condensate fraction $`N_0/N`$ versus the temperature $`T`$. In this figure we include also the quasi-classical approximation of the Popov and Hartree-Fock (i.e. to put the functions $`v_j`$ equal to zero) method and the curve of the ideal gas in the thermodynamic limit (13). The critical temperature is in good agreement with the analytical estimation obtained by using (15) and (17): the correction (15) shifts the critical temperature from $`63`$ to $`58.4`$ nK, whereas the correction (17) from $`58.4`$ to $`58.5`$ nK. It is interesting to observe that the quasi-classical approximation of the Popov method gives accurate results for the condensate fraction at all temperatures here considered. To obtain the quasi-classical curve we have followed the self-consistent procedure described in Ref. 8, which avoids the numerical diagonalization of the Bogoliubov-Popov equations by using the quasi-classical formula (10) for the excitation spectrum. It is also worth noting that in these conditions the Hartree-Fock approximation gives results compatible to those of the more accurate Popov method.

Finally, we evaluate the loss rate of the Bose condensate due to two- and three-body collisions. We use the relation :

$$\mathrm{\Gamma }_L=\frac{dN_0}{dt}=Kd^3𝐫|n_0(𝐫)|^2+Ld^3𝐫|n_0(𝐫)|^3$$
(20)   
where $`K=1.210^{14}cm^3s^1`$ is the coefficient of the two-body dipolar collisions and $`L=2.610^{28}cm^6s^1`$ is the coefficient of the three-body recombination collisions. These coefficients are strictly valid at zero temperature because they depend on temperature. Nevertheless, we can use Eq. (20) to estimate the decay rate also at finite temperature because, in our self-consistent calculations, the condensate density is modified by the presence of the thermal cloud. We expect that our results are quite accurate near $`T=0`$.

In Fig. 6 we plot, for $`N=1000`$, the mean radius of the condensate, obtained with the self-consistent Popov calculation, and the total loss rate $`\mathrm{\Gamma }_L`$ as a function of temperature. By increasing the temperature, the number $`N_0`$ of condensed atoms decreases (see Fig. 5) and, because of the attractive interaction, the radius of the condensate grows: the loss rate is reduced due to a smaller probability of inelastic collision among condensed atoms. We can proceed in a different way by keeping fixed, at all $`T`$, the number $`N_0`$ of condensed atoms and, as a consequence, by increasing the total number $`N`$ of atoms in the trap. In this case the loss rate is practically constant because there is no significant modifications in the condensate density.

### 4.3 Near the critical threshold $`N_c`$

In this section we study the behavior of the Bose vapor by changing the number of the atoms. In particular we consider a total number of atoms near the critical threshold $`N_c`$, at which there is the collapse of the condensate.

At zero temperature, by solving the Popov equations (6) and (7), we find $`N_c=1260`$, while the condensed atoms are $`N_0^c=1256`$. Remember that previously, by using the GP equation alone, i.e. neglecting the quantum depletion, we have found $`N_c=N_0^c=1287`$. This discrepancy is due to the properties of the spectrum of elementary excitations. By increasing $`N`$ there are no relevant changes in the the excitation energies besides the monopole mode $`(n=1,l=0)`$ that decreases. Strong oscillations begin when there is a level crossing between this mode and the dipole mode $`(n=0,l=1)`$ that has lower energy. At this point there are two modes with the same frequency and the system becomes unstable because of the induced resonance.

We use the quasi-classical Popov approximation to evaluate the critical threshold $`N_c`$ as a function of the temperature. We modify the numerical procedure substituting the first Bogoliubov step by a Popov step, using a Boltzmann function $`\mathrm{exp}\left(V_0(𝐫)/k_BT\right)`$ as trial non-condensed density. We have compared these quasi-classical results with the numerical ones of the full Popov approximation for some $`T`$ and $`N`$. For temperature less than about $`80`$ nK and total number of atoms less than about $`4000`$ we find a difference of about $`3\%`$ for the condensate fraction but when the temperature increases the quasi-classical method underestimates the condensed number of atoms.

In Fig. 7a we plot the condensate fraction as a function of the temperature $`T`$ for different numbers of <sup>7</sup>Li atoms. By increasing the temperature it is possible to put a larger number of atoms in the trap but the metastable condensate has a decreasing number of atoms. These results, obtained with the quasi-classical Popov approximation, show that, when $`\omega =878`$ Hz, for $`N<N_c=1260`$ the system is always metastable and it has a finite condensate fraction for $`0T<T_{_{BEC}}`$. For $`N>N_c`$ the system is metastable only for temperatures that exceed a critical value $`T_c`$ and it has a finite condensate fraction for $`T_c<T<T_{_{BEC}}`$. Note that the critical number $`N_0^c`$ of condensed atoms slightly decreases by increasing the temperature.

The collapse occurs because, by increasing the temperature, the effective potential $`V_0(𝐫)+2gn_{nc}(𝐫)`$ becomes tighter and tighter and so its characteristic length is reduced. We have seen that when the condensate fraction is only few percents then $`n_{nc}(0)/n_0(0)`$ becomes sufficiently large and the BEC cannot take place anymore. Thus, our calculations suggest that there should be a critical temperature $`T^{}`$ beyond which the BEC transition is practically inhibited, independently of the number of atoms in the trap. The data in  suggest that there is a very small fraction of condensed atoms also for temperature greater than $`400`$ nK. Actually, our rate of decrease of $`N_c^0`$ with $`T`$ is greater than the Hartree-Fock one obtained in Ref. 20. Note, however, that these mean-field methods are reliable when there is a large fraction of atoms in the condensate.

Obviously $`T_c`$ and $`T_{_{BEC}}`$ are functions of $`N`$. Such temperatures, and also the critical parameter $`N_c`$, are functions of the trap geometry, namely the frequency $`\omega `$ (see Fig. 7b), and of the scattering length $`a_s`$.

Recently, the finite-temperature behavior of condensed <sup>7</sup>Li near the critical threshold has been also studied by Davids et al . They have used a different numberical scheme to solve self-consistently the Popov equations. Their results are consistent with ours: the critical number of condensed atoms slightly decreases by increasing the temperature. In particular, their calculations show that the rate of decrease of $`N_0^c`$ with $`T`$ is greater in the Popov treatment than in the Hartree-Fock treatment used by Houbiers and Stoof .

## 5 CONCLUSIONS

We have studied the thermodynamical properties of a weakly-interacting cloud of <sup>7</sup>Li atoms confined in a harmonic trap. Due to the attractive interaction, the Bose condensed atoms remain in a metastable state until the number of atoms is lower than a critical threshold, that is about $`10^3`$ atoms at zero temperature. Far from the threshold, we have found that our Bogoliubov-Popov results give a transition temperature that is lower of about $`5`$ nK than the ideal gas value. This numerical value is in good agreement with the correcting formulas proposed in the literature . It is important to observe that, in these conditions, also the quasi-classical approximation of the Popov method gives a reliable estimation of the condensate fraction as a function of temperature. We have calculated the total loss rate (two- and three-body collisions) of the condensate: the loss rate is strongly influenced by the number of condensate atoms and by their density profile. This density profile is not affected by increased thermal fraction of atoms in the trap.

Finally, we have studied the temperature dependence of the critical threshold. We have found that the critical number of condensed atoms decreases by increasing the temperature. The condensate presents strong oscillations when there is a crossing between the dipole mode and the monopole mode, that decreases for large numbers of atoms near the threshold. Such a behavior induces the collapse of the system: in this simple model, with local attractive interaction, the wavefunction of the condensate becomes a delta function and the energy per particle goes to minus infinity. We have shown that for $`N<N_c`$ the system is always metastable and it has a finite condensate fraction for $`0T<T_{_{BEC}}`$. Instead, for $`N>N_c`$ the system is metastable only for temperatures that exceeds a critical value $`T_c`$; in this case the system has a finite condensate fraction for $`T_c<T<T_{_{BEC}}`$. Moreover, within the quasi-classical Popov approximation, our calculations suggest that there should be a temperature $`T^{}`$, beyond which the BEC transition is inhibited for any number of particles.

We believe that our calculations could be useful to select experimental conditions for which a fraction of condensed <sup>7</sup>Li atoms is detectable at thermal equilibrium.

## ACKNOWLEDGEMENTS

This work has been supported by INFM under the Research Advanced Project (PRA) on ”Bose-Einstein Condensation”.

## Figure Captions

Fig. 1. Specific quantities of condensed <sup>7</sup>Li atoms, obtained by solving the the Gross-Pitaevskii equation, as a function of the total number of atoms (N). The energies are in units of $`\mathrm{}\omega `$, where $`\omega =878`$ Hz is the frequency of the external harmonic potential. The radius is in units of the characteristic length $`a_h=(\mathrm{}/m\omega )^{1/2}=3.23`$ $`\mu `$m of the trap.

Fig. 2. a) Elementary excitations of $`10^3`$ <sup>7</sup>Li atoms at $`T=58`$ nK. Comparison of the Popov approximation with the Boboliubov one. Units as in Fig. 1. $`n`$ is the radial quantum number. b) The lowest Popov energy levels vs. temperature. $`N=10^3`$ and units as in Fig. 1.

Fig. 3. Popov calculations for $`10^3`$ <sup>7</sup>Li atoms. a) Profile of the condensate density for temperature ranging from $`0`$ to $`58`$nK. b) The same for the non-condensate density.

Fig. 4. Number of non condensed atoms ($`N_{nc}`$) for each value of the angular quantum number $`l`$, with a comparison of the Popov approximation with the Bogoliubov one. The radius is in units of the characteristic length $`a_h=(\mathrm{}/m\omega )^{1/2}=3.23`$ $`\mu `$m of the trap.

Fig. 5. Condensate fraction vs temperature for $`10^3`$ <sup>7</sup>Li atoms.

Fig. 6. Mean radius ($`r_m`$) of the condensate (top) and total loss rate (bottom) vs. the temperature for $`10^3`$ <sup>7</sup>Li atoms. Units as in Fig. 1.

Fig. 7. Condensate fraction vs temperature for different numbers of <sup>7</sup>Li atoms. Results obtained with the quasi-classical Popov approximation. Trap frequency: a) $`\omega _1=878`$ Hz; b) $`\omega _2=439`$ Hz.

| T (nK) | $`(1,0)`$ | $`(0,1)`$ | $`(10,5)`$ | $`(10,10)`$ |
| --- | --- | --- | --- | --- |
| 0 | 1.77 | 1.00 | 25.24 | 30.30 |
| 10 | 1.77 | 1.00 | 25.23 | 30.29 |
| 20 | 1.79 | 0.99 | 25.20 | 30.26 |
| 30 | 1.83 | 0.99 | 25.13 | 30.18 |
| 40 | 1.88 | 0.98 | 25.03 | 30.07 |
| 50 | 1.92 | 0.96 | 24.91 | 29.93 |
| 53 | 1.93 | 0.96 | 24.87 | 29.88 |
| 58 | 1.95 | 0.95 | 24.81 | 29.80 |
| 58.5 | 1.95 | 0.95 | 24.80 | 29.79 |

Table 1. Some excitation energies for $`N=1000`$ <sup>7</sup>Li atoms with the Popov method at different values of the temperature $`T`$. The energies, labeled by quantum numbers $`(n,l)`$, are in units of $`\mathrm{}\omega `$, where $`\omega =878`$ Hz is the frequency of the external harmonic potential.

| T (nK) | $`N_0^{(B)}`$ | $`N_0^{(P)}`$ | $`N_{nc}^{(B)}`$ | $`N_{nc}^{(P)}`$ | $`N_0^{(B)}/N`$ | $`N_0^{(P)}/N`$ |
| --- | --- | --- | --- | --- | --- | --- |
| 0 | 999.58 | 999.58 | 0.42 | 0.42 | 99.96 | 99.96 |
| 10 | 991.02 | 990.99 | 8.98 | 9.00 | 99.10 | 99.10 |
| 20 | 950.07 | 949.52 | 49.93 | 50.48 | 95.01 | 94.95 |
| 30 | 854.31 | 850.64 | 145.69 | 149.36 | 85.43 | 85.06 |
| 40 | 680.75 | 667.18 | 319.24 | 332.82 | 68.07 | 66.72 |
| 50 | 407.76 | 371.82 | 592.24 | 628.18 | 40.77 | 37.18 |
| 53 | 303.47 | 257.48 | 696.52 | 742.52 | 30.35 | 25.75 |
| 58 | 105.00 | 38.06 | 894.99 | 961.94 | 10.50 | 3.80 |
| 58.5 | 83.41 | 13.97 | 916.59 | 986.02 | 8.34 | 1.40 |

Table 2. Number of condensed ($`N_0`$) and non condensed ($`N_{nc}`$) atoms and condensate fraction $`N_0/N`$ for $`N=1000`$ <sup>7</sup>Li atoms at different values of the temperature $`T`$. <sup>(B)</sup> denotes the Bogoliubov result and <sup>(P)</sup> the Popov one.