text:
# Quantum Phase Transitions in the Shastry-Sutherland Model for 𝐒𝐫𝐂𝐮_𝟐(𝐁𝐎_𝟑)_𝟐

\[

## Abstract

We investigate the quantum phase transitions in the frustrated antiferromagnetic Heisenberg model for $`\mathrm{SrCu}_2(\mathrm{BO}_3)_2`$ by using the series expansion method. It is found that a novel spin-gap phase, which is adiabatically connected to the plaquette-singlet phase, exists between the dimer and the magnetically ordered phases known so far. When the ratio of the competing exchange couplings $`\alpha (=J^{}/J)`$ is varied, this spin-gap phase exhibits the first- (second-) order quantum phase transition to the dimer (the magnetically ordered) phase at the critical point $`\alpha _{c1}=0.677(2)`$ ($`\alpha _{c2}=0.86(1)`$). Our results shed light on some controversial arguments about the nature of the quantum phase transitions in this model.

\]

Two-dimensional (2D) antiferromagnetic quantum spin systems with the spin gap have been the subject of considerable interest. A typical compound found recently is $`\mathrm{SrCu}_2(\mathrm{BO}_3)_2`$, in which the characteristic lattice structure of the $`\mathrm{Cu}^{2+}`$ spins (see Fig. 1) stabilizes the singlet ground state. This system has been providing a variety of interesting phenomena such as the plateaus in the magnetization curve observed at $`1/3,1/4`$ and $`1/8`$ of the full moment. The spin system may be described by the 2D Heisenberg model on the square lattice with some diagonal bonds which is referred to as the Shastry-Sutherland model, as pointed out by Miyahara and Ueda. The key structure with the orthogonal dimers shown in Fig. 1 makes the system unique and particularly interesting among 2D spin-gap compounds. In this frustrated system, there may occur non-trivial quantum phase transitions when the nearest-neighbor coupling $`J`$ and the next-nearest-neighbor coupling $`J^{}`$ are varied. Albrecht and Mila discussed the possibility of a helical phase between the dimer and the magnetically ordered phases by means of the Schwinger boson mean-field theory. Recent theoretical studies, however, have suggested that there may not be such a helical phase, but the first-order phase transition occurs from the dimer to the ordered phases . Furthermore, more recent study claims that the phase transition should be of the second order with a non-trivial critical exponent $`\nu =0.45(2)`$. These controversial conclusions may come from the fact that the quantum phase transition in the Shastry-Sutherland model suffers from the strong frustration due to the competing exchange interactions $`J`$ and $`J^{}`$, and therefore a careful treatment should be necessary to figure out the correct nature of the phase transition. In particular, we have to keep in mind that such a strong frustration may possibly stabilize another spin-gap phase distinct from the dimer phase.

In this paper, by calculating the ground state energy, the staggered susceptibility and the spin gap by means of the series expansion method, we find that there should exist a novel spin-gap phase with the disordered ground state, which is stabilized by the strong frustration, between the dimer and the magnetically ordered phases. The spin-gap phase found in this paper undergoes the first- (second-) order quantum phase transition to the dimer (the ordered) phase, when the exchange couplings $`J`$ and $`J^{}`$ are varied. The existence of the new phase can resolve controversial conclusions  deduced for the quantum phase transitions in this frustrated model. We also point out that the material $`\mathrm{SrCu}_2(\mathrm{BO}_3)_2`$ lies around the phase boundary between these two spin-gap phases, which may give a natural interpretation for the 1/8-plateau formation in the magnetization curve.

To investigate the frustrated spin system for the compound $`\mathrm{SrCu}_2(\mathrm{BO}_3)_2`$, we consider the 2D quantum Heisenberg model (Shastry-Sutherland model), which is described by the following Hamiltonian

$`H=J{\displaystyle \underset{n.n.}{}}𝐒_i𝐒_j+J^{}{\displaystyle \underset{n.n.n.}{}}𝐒_i𝐒_j,`$ (1)
where $`𝐒_i`$ is the $`s=1/2`$ spin operator at the $`i`$-th site and $`J`$ $`(J^{})`$ represents the nearest-neighbor (next-nearest-neighbor) antiferromagnetic exchange coupling.

For later convenience, we introduce the ratio $`\alpha =J^{}/J`$. In Fig. 1, we have drawn the 2D Heisenberg model schematically. We note that the system with only the next-nearest-neighbor coupling $`J^{}`$ is equivalent to the Heisenberg model on the square lattice which has a spontaneous staggered magnetization at $`T=0`$. From this point of view, the nearest-neighbor coupling $`J`$ is regarded as the coupling for a diagonal bond (see Fig. 2), which gives rise to the frustration together with $`J^{}`$ .

In order to study the quantum phase transitions in this spin system, we employ the series expansion method developed by Singh, Gelfand and Huse. We recall here that the quantum phase transitions in the Shastry-Sutherland model have been discussed by Weihong et al. and Müller-Hartmann et al., by means of the dimer and the Ising expansions, from which the critical point between the dimer phase and the magnetically ordered phase has been estimated as $`\alpha _c=0.691(6)`$ and $`0.697(2)`$, respectively. As mentioned above, however, there is a controversy to be resolved about the nature of the phase transitions. Also, in order to determine the complete phase diagram, it is crucial to figure out whether there may exist another spin-gap phase besides the above two phases. We will address this problem in the following by using the series expansion method.

To see our strategy clearly, we start with the 2D quantum spin model schematically shown in Fig. 2 , which is topologically equivalent to the original model in Fig. 1.

In this figure, we have introduced an auxiliary parameter $`\lambda `$, which parameterizes the antiferromagnetic couplings labeled by the bold solid, the thin solid and the dashed lines, respectively, as $`J^{}`$, $`\lambda J^{}`$ and $`\lambda J^{}/\alpha (=\lambda J)`$. Note that the system is reduced to the original Shastry-Sutherland model in the case of $`\lambda =1`$. An important point is that the introduction of $`\lambda `$ enables us to perform the cluster expansion starting from the isolated plaquette singlets  ($`\lambda =0`$), which naturally describes the most likely spin-gap phase distinct from the dimer phase.

To proceed the analysis based on the series expansion, we divide the original Hamiltonian eq. (1) into two parts as $`H=J^{}\left[𝐒_i𝐒_j+\lambda \mathrm{\Gamma }_{ij}𝐒_i𝐒_j\right]`$, where $`\mathrm{\Gamma }_{ij}=1`$ or $`\alpha ^1`$ for each bond on the square lattice (see Fig. 2). The first term is the unperturbed Hamiltonian which stabilizes the isolated plaquette singlets with the spin excitation gap. The perturbed Hamiltonian labeled by $`\lambda `$ connects these isolated plaquette singlets, by which a 2D network develops. We expand the staggered susceptibility $`\chi _{\mathrm{AF}}`$, the spin-triplet excitation energy $`E(𝐤)`$ and the ground state energy $`E_\mathrm{g}`$ as a power series in $`\lambda `$. Here, to estimate the susceptibility, we introduce the Zeeman term $`H^{}=h\left[_{iA}S_i^z_{iB}S_i^z\right]`$, where $`h`$ is the staggered magnetic field and $`A(B)`$ denotes one of the two sublattices. Note that an asymptotic analysis of the series expansion is necessary to deduce the accurate phase boundary on which the susceptibility $`\chi _{\mathrm{AF}}`$ diverges and the spin gap $`\mathrm{\Delta }=E(𝐤=\mathrm{𝟎})`$ vanishes. For this purpose, we make use of the Padé approximants for both quantities obtained up to the finite order in $`\lambda `$. Besides ordinary Dlog Padé approximants, we also employ biased Padé approximants, for which we assume that the phase transition in our 2D quantum spin models should belong to the universality class of the 3D classical Heisenberg model. Then the critical value of $`\lambda _c`$ is determined by the formula $`\chi _{\mathrm{AF}}(\lambda _c\lambda )^\gamma `$ and $`\mathrm{\Delta }(\lambda _c\lambda )^\nu `$ with the known exponents $`\gamma =1.4`$ and $`\nu =0.71`$.

We first calculate the staggered susceptibility $`\chi _{\mathrm{AF}}`$ and the spin gap $`\mathrm{\Delta }`$ by means of the plaquette expansion up to the fourth and the fifth order in $`\lambda `$, respectively, for various values of $`\alpha `$. Using the Dlog and the biased Padé approximants, we end up with the phase diagram shown in Fig. 3.

In this figure, the solid (dashed) line represents the phase boundary obtained by the biased Padé approximants for the spin gap (the staggered susceptibility). When $`\alpha \mathrm{}`$ and $`\lambda =0`$, the system is reduced to an assembly of the isolated plaquettes with the spin gap. As $`\lambda `$ is increased, the correlation between these plaquettes grows up and the second-order quantum phase transition from the spin-gap phase to the magnetically ordered phase occurs at the critical point $`\lambda _c=0.56`$ for $`\alpha \mathrm{}`$, which has already been studied by several groups. On the other hand, decreasing $`\alpha `$ enhances the frustration, which in turn suppresses the antiferromagnetic correlation, thus shifting the phase boundary upward for smaller $`\alpha `$ in the phase diagram. It is seen that two lines obtained from the distinct quantities are in good agreement with each other, which implies that the obtained phase boundary is rather accurate in spite of the lower-order pertubative calculation. By exploiting the phase boundary determined by means of biased Padé approximants for the spin gap, the critical value is given by $`\alpha _{c2}=0.86(1)`$ for $`\lambda =1`$. Recall that the system is reduced to the original model only for $`\lambda =1`$. We thus find that the Shastry-Sutherland model has the disordered ground state in the region ($`0<\alpha <\alpha _{c2}`$) on the $`\lambda =1`$ line.

The above result does not necessarily imply that in the region $`0<\alpha <\alpha _{c2}`$ the system always belongs to the disordered phase which is continuously connected to isolated plaquettes. In fact, it is known that the orthogonal dimer phase appears in the vicinity of $`\alpha =0`$ . Therefore, it is necessary to clarify how these two spin-gap phases compete with each other by carefully comparing the ground state energy $`E_\mathrm{g}`$. To this end, performing the plaquette expansion up to the seventh order in $`\lambda `$ with $`\alpha `$ being fixed, we estimate the ground state energy $`E_\mathrm{g}`$ for the Shastry-Sutherland model $`(\lambda =1)`$ by means of the first order inhomogeneous differential method. The results are shown in Fig. 4.

As mentioned above, the system stabilizes the orthogonal dimer ground state for smaller $`\alpha `$. It is found, however, that the first-order transition to the novel spin-gap phase introduced here occurs at the critical point $`\alpha _{c1}=0.677(2)`$. It is also seen from this figure that further increase of $`\alpha `$ induces the antiferromagnetic order, whose transition point is determined by the crossing point of the ground-state energy obtained respectively by the Ising and plaquette expansions. The result confirms the second-order phase transition deduced above, and the transition point estimated from the figure is consistent with $`\alpha _{c2}=0.86(1)`$ obtained by the analysis of the susceptibility and the spin gap. Consequently, we end up with the phase diagram for the Shastry-Sutherland model as shown in Fig. 5. The present results shed light on the controversial arguments whether the quantum phase transition in this model is of the first or second order. In those previous studies, it was believed that the phase transition occurs only once between the dimer phase (I) and the ordered phase (III), giving rise to some confusions. Our phase diagram clearly resolves this problem by explicitly showing the existence of the new spin-gap phase (II) which undergoes the first- (I$``$II) as well as the second-order transitions (II$``$III).

To check the validity of the above phase diagram, we also show the results for the spin gap as a function of $`\alpha =J^{}/J`$ in Fig. 6.

In this figure, the results obtained by Weihong et al. are shown for the orthogonal dimer phase (I: $`0<\alpha <\alpha _{c1}`$). In the new phase (II: $`\alpha _{c1}<\alpha <\alpha _{c2}`$), we determine the values of the spin gap at $`𝐤=\mathrm{𝟎}`$ by means of the plaquette expansion up to the fifth order in $`\lambda `$ with the first order inhomogeneous differential method. The results are shown as the dots with the error bars. As seen in this figure, with the decrease of $`\alpha `$ from the second-order transition point $`\alpha _{c2}`$, the spin-gap continuously grows up to stabilize the disordered ground state. As $`\alpha `$ is further decreased, the first-order phase transition occurs at $`\alpha _{c1}`$.

In order to further confirm the present results, we have performed a different series expansion by choosing the isolated plaquettes with diagonal bonds as an initial configuration, which is different from the one shown in Fig. 2. The calculation of the susceptibility up to the fourth order yields second-order transition with $`\alpha _{c2}=0.87(3)`$, being consistent with the above results. Furthermore, to confirm the first-order phase transition between the two spin-gap states, we have checked how the first-order phase transition point known for the 1D orthogonal-dimer chain  evolves with the increase of the inter-chain couplings. By performing the exact diagonalization studies for the $`4\times 4`$ system, we have found that the first-order transition point for 1D is continuously changed, and in the Shastry-Sutherland case, it coincides with the one found above within reasonable accuracy $`(\alpha _{c1}0.66)`$, providing further support to our conclusion on the phase diagram. Although our results still seem to be partly contradicted to the staggered magnetization obtained by Weihong et al., we believe that this could be resolved by further analysis of the results of the Ising expansion.

Before concluding the paper, a brief comment is in order for the plateau-formation in the magnetization curve. Experimentally, the plateaus in the magnetization curve have been observed for the compound $`\mathrm{SrCu}_2(\mathrm{BO}_3)_2`$ at 1/3, 1/4 and $`1/8`$ of the full moment. In the theoretical studies on the dimer phase, it has been clarified that the stripe order of the isolated dimer-triplets is important to understand the $`1/3`$ and $`1/4`$plateaus. On the other hand, it is not so trivial why the 1/8-plateau occurs in this compound, although a possible mechanism has been proposed . We think that the formation of the 1/8 plateau may reflect the fact that this compound is located around the first-order phase transition point between the two spin-gap phases and thereby possesses the dual properties inherent in two distinct phases implicitly. We note here that the new spin-gap phase belongs to the same phase as the Heisenberg model on the 1/5-depleted square lattice proposed for $`\mathrm{CaV}_4\mathrm{O}_9`$. Therefore it is likely that the 1/8-plateau could occur in the same origin discussed by Momoi and Totsuka for the plaquette system related to the 1/5-depleted Heisenberg model. It is interesting to further clarify the mechanism of the 1/8-plateau by taking into account the above dual properties explicitly, which is now under consideration.

In conclusion we have discussed the phase diagram for the Shastry-Sutherland model for the compound $`\mathrm{SrCu}_2(\mathrm{BO}_3)_2`$ by means of the series expansion method. Our analysis has shown that there exists a novel spin-gap phase with the disordered ground state, which is adiabatically connected to the plaquette-singlet phase, between the dimer and the magnetically ordered phases known so far. When the exchange coupling ratio $`\alpha =J^{}/J`$ is varied, the first-order phase transition occurs from the dimer state to the new spin-gap state, while the second-order phase transition occurs from this spin-gap state to the magnetically ordered state. This sheds light on the nature of the quantum phase transitions in this model, and resolves apparently controversial conclusions on this issue.

We would like to thank K. Ueda, S. Miyahara and K. Okunishi for useful discussions. The work is partly supported by a Grant-in-Aid from the Ministry of Education, Science, Sports, and Culture. A. K. is supported by the Japan Society for the Promotion of Science. A part of numerical computations in this work was carried out at the Yukawa Institute Computer Facility.


format: text

text:
# Kinetic Theory of Collective Excitations and Damping in Bose-Einstein Condensed Gases

## I Introduction

The nonzero-temperature dependence of the frequencies and damping of collective modes in trapped atomic Bose gases has been investigated extensively both experimentally  and theoretically . Above the Bose-Einstein transition temperature the lowest-lying collective modes have been calculated in the collisionless regime using a collisionless Boltzmann or Landau-Vlasov equation , and in the hydrodynamic regime either by using the conservation laws of hydrodynamics  or by using the appropriate quantum kinetic equation . In the collisionless regime the frequency of the collective mode is large compared to the mean collision frequency, in contrast to the situation in the hydrodynamic regime. Several papers have also studied the regime intermediate between the collisionless and hydrodynamic regimes by taking into account interatomic collisions , which turn out to mostly lead to mode damping. These papers have indicated that the experiments were performed under conditions intermediate between collisionless and hydrodynamic. (This is particularly true for the experiments performed by Mewes et al.  and Stamper-Kurn et al. ). For temperatures far below the transition temperature the collisionless modes can be described accurately by the time-dependent Gross-Pitaevskii equation . At higher temperatures the noncondensate fraction becomes substantial and the modes of the condensate are now coupled to those of the thermal cloud. In Refs.  the temperature dependence of the mode frequencies has been calculated by employing the Popov approximation to include the static mean-field effects of the noncondensate atoms in the Gross-Pitaevskii equation. This was improved later by taking into account also the dynamics of the thermal cloud using a RPA approximation  or the collisionless Boltzmann equation . For the hydrodynamic regime a two-fluid model has been developed in Refs. . Moreover, the theory of Zaremba, Griffin, and Nikuni  was improved later by the same authors to include collisions between condensate and noncondensate atoms .

Using a combination of two previous papers, namely Refs.  and , we aim in this paper at interpolating between the collisionless and hydrodynamic regimes for experiments below the transition temperature. The authors of Ref.  have already performed such an interpolation for a gas above the critical temperature by using a Boltzmann equation with a relaxation time approximation. In contrast, the authors of Ref.  use a collisionless Boltzmann equation for the thermal cloud coupled to a time-dependent nonlinear Schrödinger equation for the condensate to consider the collisionless dynamics below the critical temperature. Here, we thus combine these two approaches by using again the Boltzmann equation and time-dependent nonlinear Schrödinger equation, but now including also the effects of interatomic collisions in the manner as put forward in Refs.  and . This implies that we have to add two collision terms to the Boltzmann equation. The first collision term represents collisions between two noncondensate atoms and the second describes collisions between a condensate and a noncondensate atom. We use for both these collision terms a relaxation time approximation, since this approximation leads above the transition temperature to a good agreement with microscopic calculations as well as with experimental data . Furthermore, for consistency reasons we also have to include a damping term in the time-dependent non-linear Schrödinger equation, which is due to collisions between the condensate and noncondensate atoms. As a result our calculation will essentially be characterized by two parameters, namely $`\tau _{22}`$ and $`\tau _{12}`$, denoting the mean collision time for collisions between the noncondensate atoms and between condensate and noncondensate atoms, respectively. Note that we use here the same notation as in Ref. . This will allow us to investigate the collisionless and hydrodynamic limits with respect to both $`\tau _{22}`$ and $`\tau _{12}`$, and enable us to fit our results for frequencies and damping with the experimental data.

To solve the complicated nonlinear dynamics of the gas, we employ a gaussian trial function for the wave function of the condensate atoms with three complex time-dependent variational parameters. For the thermal cloud we use a distribution function that incorporates deviations from the Bose-Einstein distribution function. The deviation function is a truncated power expansion in the momenta and coordinates of the atoms. Our system of coupled equations is then assembled from the Euler-Lagrange equations for the variational parameters in the condensate wavefunction and from all the linear and quadratic moments of the Boltzmann equation for the distribution function of the noncondensate atoms. The solution of the linearized version of this system of coupled equations results in a dispersion relation that gives the frequencies and damping of the coupled modes in terms of $`\tau _{22}`$ and $`\tau _{12}`$. Although we restrict our calculation to the axially symmetric traps used in the experiments, the generalization to fully anisotropic traps is straightforward.

The rest of this paper is organized as follows: In the next section we present the theoretical details of our calculation. Specifically, we write down the time-dependent nonlinear Schrödinger equation for the condensate wave function and the linearized Boltzmann equation for the thermal cloud, and show how to treat the two collision terms in a relaxation time approximation. In section II A we present our trial functions for the condensate wave function and noncondensate distribution function and in section II B we discuss how we obtain the equilibrium state of the gas around which we have to expand to find the collective modes. In section III we present the resulting dispersion relation and discuss its collisionless and hydrodynamic limits. We then compare our results with experiment. We end in section IV by summing up our main conclusions.

## II Coupled Dynamics of the condensate and the thermal cloud

It has been shown that the coupled time-dependent nonlinear Schrödinger equation and the quantum Boltzmann equation can both be derived starting from the equation of motion

$`i\mathrm{}{\displaystyle \frac{\widehat{\psi }(𝐫,t)}{t}}`$ $`=`$ $`\text{(}{\displaystyle \frac{\mathrm{}^2^2}{2m}}+V^{ext}(𝐫)`$ (1)
$`+`$ $`g\widehat{\psi }^{}(𝐫,t)\widehat{\psi }(𝐫,t)\text{)}\widehat{\psi }(𝐫,t),`$ (2)
for the Heisenberg field operator $`\widehat{\psi }(𝐫,t)`$. Here, $`g`$ is the effective two-body interaction, which is given in terms of the scattering length $`a`$ and the atomic mass $`m`$ as $`g=4\pi a\mathrm{}^2/m`$. The external potential we take here is a harmonic potential of the general form,

$$V^{ext}(𝐫)=\frac{1}{2}m(\omega _1^2x^2+\omega _2^2y^2+\omega _3^2z^2),$$
(3)
where $`\omega _i`$ is the characteristic frequency of the trap in the $`i`$th direction. To obtain the time-dependent nonlinear Schrödinger equation we write $`\widehat{\psi }(𝐫,t)`$ as

$$\widehat{\psi }(𝐫,t)=\mathrm{\Phi }(𝐫,t)+\widehat{\psi }^{}(𝐫,t),$$
(4)
where $`\mathrm{\Phi }(𝐫,t)`$ is the appropriate nonequilibrium expectation value of $`\widehat{\psi }(𝐫,t)`$ and the operator $`\widehat{\psi }^{}(𝐫,t)`$ describes the noncondensate atoms. The number density of condensate atoms $`n_c(𝐫,t)`$ is related to $`\mathrm{\Phi }(𝐫,t)`$ by $`n_c(𝐫,t)=|\mathrm{\Phi }(𝐫,t)|^2`$ while the number density of the noncondensate atoms $`n^{}`$ equals $`\widehat{\psi ^{}}^{}(𝐫,t)\widehat{\psi ^{}}(𝐫,t)`$, where the brackets represent the averaging over the nonequilibrium density matrix. Substituting this form for $`\widehat{\psi }(𝐫,t)`$ in the equation of motion and averaging, we obtain

$`i\mathrm{}{\displaystyle \frac{\mathrm{\Phi }(𝐫,t)}{t}}`$ $`=`$ $`\text{(}{\displaystyle \frac{\mathrm{}^2^2}{2m}}+V^{ext}(𝐫)+gn_c(𝐫,t)+2gn^{}(𝐫,t)`$ (5)
$``$ $`iR(𝐫,t){\displaystyle \frac{}{}}\text{)}\mathrm{\Phi }(𝐫,t),`$ (6)
where $`R(𝐫,t)`$ is given by

$`R(𝐫,t)`$ $`=`$ $`g(\widehat{\psi ^{}}(𝐫,t)\widehat{\psi ^{}}(𝐫,t)\mathrm{\Phi }^{}(𝐫,t)`$ (7)
$`+`$ $`\widehat{\psi ^{}}^{}(𝐫,t)\widehat{\psi ^{}}(𝐫,t)\widehat{\psi ^{}}(𝐫,t))/\mathrm{\Phi }(𝐫,t).`$ (8)
Next, the quantum Boltzmann equation can be derived by writing an equation of motion for the distribution function of the noncondensate atoms $`f(𝐩,𝐫,t)`$. This is usually done by writing $`f(𝐩,𝐫,t)`$ as a Wigner transform of $`\widehat{\psi ^{}}^{}(𝐫^{},t)\widehat{\psi ^{}}(𝐫,t)`$ and then determining the time evolution of $`f(𝐩,𝐫,t)`$ from the equation of motion in Eq. (2) . In the Hartree-Fock approximation the resulting Boltzmann equation takes the form

$$\left[\frac{f}{t}+_𝐩E_𝐫_𝐫E_𝐩\right]f=C_{22}[f]+C_{12}[f],$$
(9)
where $`C_{22}`$ is the contribution to the rate of change of $`f`$ due to collisions between noncondensate atoms, while $`C_{12}`$ is the contribution due to collisions between the condensate and the noncondensate atoms. The energy $`E(𝐩,𝐫,t)`$ of the noncondensate atoms is in this approximation given by

$$E(𝐩,𝐫,t)=p^2/2m+V^{ext}(𝐫)+2gn(𝐫,t),$$
(10)
where $`n=n_c+n^{}`$ is the total density.

In Ref.  explicit forms of the two collision terms have been written down as follows

$`C_{22}[f]`$ $`=`$ $`{\displaystyle \frac{2g^2}{(2\pi )^5\mathrm{}^7}}{\displaystyle 𝑑𝐩_2𝑑𝐩_3𝑑𝐩_4\delta (𝐩+𝐩_2𝐩_3𝐩_4)}`$ (11)
$`\times `$ $`\delta (E+E_2E_3E_4)`$ (12)
$`\times `$ $`\left[(1+f)(1+f_2)f_3f_4ff_2(1+f_3)(1+f_4)\right],`$ (13)
and

$`C_{12}[f]`$ $`=`$ $`{\displaystyle \frac{2g^2n_c}{(2\pi )^2\mathrm{}^4}}{\displaystyle 𝑑𝐩_1𝑑𝐩_2𝑑𝐩_3}`$ (14)
$`\times `$ $`\delta (m𝐯_c+𝐩_1𝐩_2𝐩_3)`$ (15)
$`\times `$ $`\delta (E_c+E_1E_2E_3)`$ (16)
$`\times `$ $`\left[\delta (𝐩𝐩_1)\delta (𝐩𝐩_2)\delta (𝐩𝐩_3)\right]`$ (17)
$`\times `$ $`\left[(1+f_1)f_2f_3f_1(1+f_2)(1+f_3)\right].`$ (18)
Here $`E_c(𝐫,t)`$ and $`𝐯_c(𝐫,t)`$ are the local energy and velocity of the condensate atoms, $`E_i`$ is a shorthand for $`E(𝐩_i,𝐫,t)`$, and similarly $`f_i`$ is a shorthand for $`f(𝐩_i,𝐫,t)`$. As is shown in Ref. , the conservation of the total number of atoms constraints the collision term $`C_{12}`$ in the Boltzmann equation to be related to the damping term $`R(𝐫,t)`$ in the complex nonlinear Schrödinger equation as

$$2n_c(𝐫,t)R(𝐫,t)=\frac{d𝐩}{(2\pi )^3}C_{12}[f].$$
(19)
The coupled equations given in Eqs. (6), (9), and (19) in principle fully determine the dynamics of the Bose-Einstein condensed gas in the Hartree-Fock approximation. Since the dynamics of the thermal cloud is experimentally only important for temperatures which are larger or comparable to the mean-field interactions we believe that this is sufficiently accurate and we do not need to use the Popov approximation at first instance. Nevertheless, these equations are too complicated to be solved exactly and some approximation is called for.

### A Trial Functions

For the noncondensate atoms we start by linearizing the Boltzmann equation in small deviations of the distribution function around its equilibrium value, namely

$$f(𝐩,𝐫,t)=f^{(0)}(𝐩,𝐫)+f^{(0)}(𝐩,𝐫)(1+f^{(0)}(𝐩,𝐫))\psi (𝐩,𝐫,t),$$
(20)
where $`f^{(0)}(𝐩,𝐫)`$ is the equilibrium distribution function given by

$$f^{(0)}(𝐩,𝐫)=\left\{\mathrm{exp}\left[(E(𝐩,𝐫)\mu )/k_BT\right]1\right\}^1.$$
(21)
Here $`\mu `$ is the chemical potential and $`k_B`$ is Boltzmann’s constant. As a result the linearized Boltzmann equation takes the form

$$\left[\frac{}{t}+_𝐩E_𝐫_𝐫E_𝐩\right]\psi =C_{22}[\psi ]+C_{12}[\psi ].$$
(22)
It is shown in Ref.  that for an uncondensed gas a trial function for $`\psi (𝐩,𝐫,t)`$ of the form

$`\psi `$ $`=`$ $`A_1x^2+B_1xp_x+C_1p_x^2+A_2y^2+B_2yp_y`$ (23)
$`+`$ $`C_2p_y^2+A_3z^2+B_3zp_z+C_3p_z^2,`$ (24)
where $`A_i(t)`$, $`B_i(t)`$, and $`C_i(t)`$ are nine time-dependent functions, is appropriate to describe the low-lying collisionless breathing modes with frequencies $`2\omega _i`$, some of which have been observed experimentally. We will therefore assume the above expansion to be also reasonably accurate below the critical temperature. For the collision integrals $`C_{22}`$ and $`C_{12}`$ we use a relaxation time approximation. In such an approximation one associates these collision integrals to mean relaxation times. However, the approximate expressions for $`C_{22}`$ and $`C_{12}`$ should still take into account the conservation laws associated with the collision processes exactly. In the case of $`C_{22}`$, which represents collisions between the noncondensate atoms, the number of atoms in the thermal cloud, their total momentum and their total kinetic energy should all be conserved. In our ansatz for $`\psi `$, terms like $`x^2`$ and $`xp_x`$ should therefore not be affected by such collisions since they correspond to two collision invariants, namely the number of atoms and the total momentum in the $`x`$ direction, respectively. On the other hand, terms like $`p_x^2`$ will be affected by collisions since it is only the sum $`p^2=p_x^2+p_y^2+p_z^2`$ which is conserved during the collision process. Therefore, we can write the following expression for the linear operator $`C_{22}`$

$`C_{22}[\psi _i]`$ $`=`$ $`{\displaystyle \frac{1}{\tau _{22}}}\{\begin{array}{ccc}(\psi _ip^2/3)& ,& \psi _i=p_x^2,p_y^2,p_z^2\\ 0& ,& \mathrm{otherwise}\end{array},`$ (27)
where $`\tau _{22}`$ is a mean collision time for the noncondensate-noncondensate atomic collisions. Note that in this expression $`C_{22}[p^2]=0`$, ensuring the conservation of the total kinetic energy.

For $`C_{12}`$, which represents collisions between the condensate and the noncondensate atoms, the number of atoms is not conserved, since the collision process involves transport of atoms back and forth from the condensate into the thermal cloud. This statement means mathematically that the zeroth moment of $`C_{12}`$, i.e., $`𝑑𝐩C_{12}`$, does not vanish in contrast to the case of $`C_{22}`$ where $`𝑑𝐩C_{22}=0`$. As a first attempt to associate with $`C_{12}`$ a mean collision time $`\tau _{12}`$ for noncondensate-condensate atomic collisions, we may write $`C_{12}1/\tau _{12}`$. However, we observe from Eq. (18) that $`C_{12}n_c(𝐫)`$. This dependence on $`n_c(𝐫)`$ requires us to assign to $`\tau _{12}`$ a position dependence that follows from $`n_c(𝐫)`$ as

$$\frac{1}{\tau _{12}(𝐫)}=\frac{1}{\tau _{12}(\mathrm{𝟎})}\frac{n_c(𝐫)}{n_c(\mathrm{𝟎})},$$
(29)
where $`\tau _{12}(\mathrm{𝟎})`$ and $`n_c(\mathrm{𝟎})`$ are the noncondensate-condensate mean atomic collision time and condensate density at the center of the trap, respectively. Generally, $`\tau _{12}(\mathrm{𝟎})`$ may, just like $`\tau _{22}`$, also have temperature dependence. By multiplying and dividing $`1/\tau _{12}(𝐫)`$ by $`\mathrm{}g`$ we can thus take $`C_{12}`$ as obeying

$$C_{12}\alpha \frac{2gn_c(𝐫)}{\mathrm{}},$$
(30)
where we introduced the dimensionless constant

$$\alpha =\frac{\mathrm{}}{2gn_c(\mathrm{𝟎})\tau (\mathrm{𝟎})}.$$
(31)
This form for $`C_{12}`$ is convenient for a reason that becomes clear lateron in this section. It should be emphasized here that the dimensionless parameter $`\alpha `$ is the second free parameter in our phenomenological calculation. The first dimensionless parameter being $`1/\overline{\omega }\tau _{22}`$ where $`\overline{\omega }=\sqrt[3]{\omega _1\omega _2\omega _3}`$. Therefore, the collisionless and hydrodynamic regimes, as well as the intermediate regime, will be described by only these two parameters. Finally, we make the relation between $`C_{12}`$ and $`\alpha `$ precise by taking

$`C_{12}[\psi _i]`$ $`=`$ $`\alpha {\displaystyle \frac{2gn_c(𝐫)}{\mathrm{}}}`$ (32)
$`\times `$ $`\{\begin{array}{ccc}(\psi _ip^2/3)& ,& \psi _i=p_{x}^{}{}_{}{}^{2},p_{y}^{}{}_{}{}^{2},p_{z}^{}{}_{}{}^{2}\\ \psi _i& ,& \psi _i=x^2,y^2,z^2\\ 0& ,& \mathrm{otherwise}\end{array}.`$ (36)
This completes our description of the treatment of the thermal cloud. Next we have to consider the condensate.

It is known that a time-dependent gaussian ansatz for the condensate wave function gives the correct frequencies of the lowest modes at zero temperature . Furthermore, the gaussian ansatz has also been used in Ref.  for the whole temperature range below the transition temperature and leads to rather good agreement with experiment. Therefore, we employ here again a gaussian ansatz for the wave function of the condensate. It has the following form

$`\mathrm{\Phi }(𝐫,t)`$ $`=`$ $`\sqrt[4]{{\displaystyle \frac{8N_cb_{1r}b_{2r}b_{3r}}{\pi ^3}}}\mathrm{exp}\left[(b_1x^2+b_2y^2+b_3z^2)\right],`$ (37)
where $`b_1`$, $`b_2`$, and $`b_3`$ are complex time-dependent variational parameters and $`b_{1r}`$, $`b_{2r}`$, and $`b_{3r}`$ are their real parts, respectively. Similarly, we denote the imaginary parts, which will appear lateron, as $`b_{1i}`$, $`b_{2i}`$, and $`b_{3i}`$. The prefactor of $`\mathrm{\Phi }(𝐫,t)`$ guarantees its normalization $`𝑑𝐫|\mathrm{\Phi }(𝐫,t)|^2`$ to be equal to the number of condensate atoms $`N_c`$.

To obtain the first set of our coupled equations of motion, we start by writing down the energy functional that corresponds to the nonlinear Schrödinger equation in Eq. (6), namely

$`E_c[\mathrm{\Phi }]`$ $`=`$ $`{\displaystyle 𝑑𝐫\text{[}\frac{\mathrm{}^2\left|𝚽\right|^2}{2m}}+V^{ext}n_c+{\displaystyle \frac{1}{2}}g(n_{c}^{}{}_{}{}^{2}+4n_cn^{})`$ (39)
$``$ $`iRn_c{\displaystyle \frac{}{}}\text{]}`$ (40)
with $`n_c=|\mathrm{\Phi }|^2`$. We evaluate this energy functional using the expression for $`\mathrm{\Phi }`$ from Eq. (LABEL:phi), and by noting that the density of the thermal cloud

$$n^{}=\frac{d𝐩}{(2\pi )^3}f(𝐩,𝐫,t),$$
(41)
and the damping term $`R(𝐫,t)`$ can be calculated using Eqs. (19) and (36). We notice here that $`iR(𝐫,t)n_c`$ turns out to be proportional to the Hartree-Fock interaction term in Eq. (40), with a proportionality constant equals to $`i\alpha `$. This explains the reason for writing $`C_{12}`$ as in Eq. (30). The equations of motion for the condensate dynamics are now the Euler-Lagrange equations resulting from varying the lagrangian

$$L=\frac{1}{2}i\mathrm{}𝑑𝐫\left(\mathrm{\Phi }^{}\frac{}{t}\mathrm{\Phi }\mathrm{\Phi }\frac{}{t}\mathrm{\Phi }^{}\right)E_c[\mathrm{\Phi }]$$
(42)
with respect to the 6 variational parameters $`b_{kr}`$ and $`b_{ki}`$.

The second set involves the equations of motion for the constants $`A_i`$, $`B_i`$, and $`C_i`$. It is obtained from taking appropriate moments of the Boltzmann equation in Eq. (22). In detail, the moments are calculated by multiplying this equation by $`f^{(0)}(𝐩,𝐫)`$ and the various components of $`\psi `$ in Eq. (24), and then integrating over $`𝐩`$ and $`𝐫`$. This results in 9 equations of motion. In combination with the previous 6 ones, we thus have 15 coupled equations of motion. The coupling is provided on the one hand by the mean-field interaction $`2gn_cn^{}`$ and the imaginary damping term $`iRn_c`$ in the condensate energy functional $`E_c`$ in Eq. (40), and on the other hand by the contribution $`2gn_c`$ to the Hartree-Fock energy $`E`$ in Eq. (10) and the $`C_{12}`$ collision term in the Boltzmann equation.

### B Equilibrium

Finally, it is important to note that linearization of the equations of motion of the condensate is obligatory to be consistent with the equations of motion of the noncondensate part which are already linearized. Therefore, we need to calculate the equilibrium state of the gas. In principle, to obtain equilibrium properties we should minimize the free energy $`F=E_{tot}TS`$, where $`E_{tot}`$ is the total energy and $`S`$ is the entropy of the gas, with respect to some variational parameters that characterize the widths of the condensate and the thermal cloud. However, a simplified estimate of the total energy and the entropy contribution to the free energy above the transition temperature shows that for the experimental conditions of interest the former is dominant over the later. This simplified estimate can be made by using a one-parameter gaussian ansatz for the thermal cloud density, namely $`n^{}\mathrm{exp}(x^2/R^2)`$, where $`R`$ is the radius of the cloud. It turns out that in the Thomas-Fermi limit the entropy contribution to the free energy is of order $`k_BT\mathrm{ln}\gamma ^{3/5}`$, where $`\gamma =N^{}a/\overline{a}`$, $`N^{}`$ is the total number of atoms in the thermal cloud, and $`\overline{a}=\sqrt{\mathrm{}/m\overline{\omega }}`$ is the harmonic oscillator length, while the total energy contribution is of order $`k_BT\gamma ^{2/5}`$. Therefore, if $`\gamma 1`$ the entropy contribution is much less than the total energy contribution.

As we shall see in the next section, the experiments were performed with a temperature-dependent total number of atoms ranging from about $`6000`$ atoms at zero temperature to approximately $`40000`$ atoms at the transition temperature. If we extend this simple estimate to temperatures below the transition temperature it turns out that the condition $`\gamma 1`$ is satisfied only for high temperatures but not for low temperatures since the number of thermal atoms becomes small. However, this condition is then no longer important since the free energy of the thermal cloud at such temperatures is small compared to the condensate energy. Therefore, we will in the following only minimize the total energy $`E_{tot}`$ with respect to the variational parameters of the condensate and the thermal cloud. To be able to calculate $`E_{tot}`$ we assume that the distribution function of the noncondensate atoms has the same form of that of a noninteracting gas but with varying spatial widths that effectively take into account the mean-field effects of both the noncondensate and condensate atoms. This effective distribution function is written as

$`f_{eff}^{(0)}`$ $`=`$ $`\left[\mathrm{exp}\left({\displaystyle \frac{p^2/2m\mu }{k_BT}}+{\displaystyle \frac{x^2}{R_{1}^{}{}_{}{}^{2}}}+{\displaystyle \frac{y^2}{R_{2}^{}{}_{}{}^{2}}}+{\displaystyle \frac{z^2}{R_{3}^{}{}_{}{}^{2}}}\right)1\right]^1,`$ (43)
where $`R_1`$, $`R_2`$, and $`R_3`$ are the widths of the noncondensate cloud in the three directions. The total energy is now a function of 6 variational parameters:

$$E_{tot}=E_{tot}(b_{1r}^{(0)},b_{2r}^{(0)},b_{3r}^{(0)},R_1,R_2,R_3),$$
(45)
where $`b_{1r}^{(0)}`$, $`b_{2r}^{(0)}`$, and $`b_{3r}^{(0)}`$ are the equilibrium values of $`b_{1r}`$, $`b_{2r}`$, and $`b_{3r}`$, respectively. The equilibrium is obtained by minimizing this energy with respect to these variational parameters.

The results of such a minimization will be shown in the next section, where we present the dispersion relation that results from solving the above-described system of 15 linearized coupled equations of motion. We discuss also the collisionless and hydrodynamic limits of these results.

## III The Dispersion Relation: Frequencies and Damping Rates

Our calculation accounts, in fully anisotropic traps, for 9 modes of the gas. In axially symmetric traps this number reduces to 6 modes. These are the in-phase and out-of-phase combinations of the two monopole ($`m=0`$) modes and one quadrupole ($`m=2`$) mode of both the condensate and the thermal cloud. Here we denote with $`m`$ the projection of the angular momentum of the mode along the axis of symmetry of the trap . We focus in this paper on the two lowest-lying $`m=0`$ and $`m=2`$ modes observed experimentally. It turns out that for the experimentally relevant temperature range, the in-phase $`m=0`$, and $`m=2`$ modes correspond mostly to oscillations of the thermal cloud, whereas the out-of-phase modes are mostly condensate oscillations. Therefore, we shall often refer in this paper to the in-phase modes as the thermal cloud or noncondensate modes and the out-of-phase modes as the condensate modes. Although our calculation provides results for another, higher-lying, monopole mode, we shall not discuss it further here. Moreover, we emphasize that throughout the following we perform our calculations for parameters taken from the experiments of Jin et al. , i.e., with <sup>87</sup>Rb atoms in an axially-symmetric trap with anisotropy ratio $`\omega _3/\omega _1=\omega _3/\omega _2=\sqrt{8}`$.

In that particular experiment the measurements were performed with a temperature-dependent total number of atoms as a result of the loss of atoms during evaporative cooling. The total number of atoms $`N_{tot}`$, as well as the number of condensate atoms $`N_c`$, are measured in the temperature range $`T/T_{BEC}0.48`$ to $`T/T_{BEC}1.0`$, where $`T_{BEC}`$ is the Bose-Einstein transition temperature. These measurements can be easily fitted with polynomials in $`T/T_{BEC}`$ as shown by the dashed line in Fig. 1. An extrapolation of such a fit to temperatures below $`0.48T_{BEC}`$ leads, however, to nonphysical situations. In particular, the two curves for $`N_{tot}`$ and $`N_c`$ cross at least once before reaching zero temperature. We overcome this problem by fitting only $`N_{tot}`$ with the experimental data, and then using $`N_c`$=$`N_{tot}(1(T/T_{BEC})^2)`$ in analogy with the ideal gas relation $`N_c`$=$`N_{tot}(1(T/T_{BEC})^3)`$ . Fig. 1 shows with the solid line the results of this slightly less accurate fit. It should be noted that using the ideal gas relation will grossly overestimate all the experimental points of $`N_c`$ for $`T/T_{BEC}>0.7`$. It turns out that the difference between the calculated frequencies using the less accurate fit from those calculated using the best fit are much smaller than the uncertainties in the measured frequencies for the experimental range of temperatures. This is shown in Fig. 2. Since we want to show also results for the complete temperature interval from zero to $`T_{BEC}`$, we employ from now on always the former fit, i.e., the solid lines in Fig. 1.

Before starting with calculating the frequencies and damping rates of the collective modes, we show in Fig. 3 the result of the minimization of the total energy that is required to obtain the equilibrium conditions of the gas, as described in the previous section. We plot the equilibrium widths of both condensate and noncondensate clouds as a function of temperature. We notice from this figure that at zero temperature and at the transition temperature we obtain rather good estimates for the widths of the condensate and thermal clouds. At zero temperature, where the thermal cloud is absent, we observe that the radial width of the condensate equals approximately $`2a_1`$, a result that is consistent with the familiar property that due to the mean-field interactions the equilibrium condensate width is larger than the ideal gas result $`a_1`$. The axial width, which is roughly $`1/\lambda ^{1/2}`$ of the radial width, is suppressed due to the anisotropy of the trap. We can also see in this figure that at the transition temperature the radial width of the condensate is slightly larger than the ideal gas one. This slight expansion is caused by the presence of the thermal cloud, whose mean-field interaction has the effect of slightly reducing the spring constants of the effective trapping potential.

Returning to the problem of the collective mode frequencies and damping, we start by neglecting collisions between the condensate and noncondensate atoms, which means that $`\alpha =0`$. In this case the dispersion relation turns out to have the following general structure

$$iw\left(P_C^{(0)}(\omega \tau _{22})^2+iP_I^{(0)}\omega \tau _{22}P_H^{(0)}\right)=0,$$
(46)
where $`P_C^{(0)}`$, $`P_I^{(0)}`$, and $`P_H^{(0)}`$ are 6th order polynomials in $`\omega ^2`$ which can all be factorized as

$$P_C^{(0)}=\underset{k=1}{\overset{6}{}}(\omega ^2\omega _{Ck}^2),$$
(47)
$$P_I^{(0)}=\underset{k=1}{\overset{6}{}}(\omega ^2\omega _{Ik}^2),$$
(48)
and

$$P_H^{(0)}=\underset{k=1}{\overset{6}{}}(\omega ^2\omega _{Hk}^2).$$
(49)
Here, $`\omega _{Ck}`$, $`\omega _{Hk}`$, and $`\omega _{Ik}`$, $`k=1,2,3`$, are temperature-dependent collisionless, intermediate, and hydrodynamic frequencies, respectively. The superscripts in $`P_C^{(0)}`$, $`P_I^{(0)}`$, and $`P_H^{(0)}`$ indicate the value of $`\alpha `$.

The general structure of this dispersion relation agrees with the result of Refs.  above the transition temperature, where this dispersion relation was studied in detail. Particularly, it was shown that damping rates calculated using this relation agree in order of magnitude with experiments and numerical calculations.

The collisionless regime is defined by the condition $`\omega \tau _{22}1`$. Therefore, we find from Eq. (46) that the collisionless frequencies are $`\omega _{Ck}`$. In the hydrodynamic limit $`\omega \tau _{22}1`$, and the hydrodynamic frequencies are $`\omega _{Hk}`$. Note that we use the word hydrodynamic here to denote that the thermal cloud is in the hydrodynamic regime though $`\alpha =0`$ and there are therefore no collisions between condensate and noncondensate atoms. In a sense this regime is thus precisely the limit discussed by Nikuni, Zaremba, and Griffin . Using the experimental parameters we calculate these frequencies for the whole temperature range below $`T_{BEC}`$. In Fig. 4 we present the results of this calculation for the $`m=0`$ and $`m=2`$ modes together with the experimental data. In this figure the collisionless curves agree with those of Bijlsma and Stoof  for $`T/T_{BEC}>0.2`$, which is not surprising since we use similar ansatz functions. The discrepancy for $`T/T_{BEC}<0.2`$ is due to the different ways of treating the equilibrium state.

Above the transition temperature, the explanation of the mode damping was based on the fact that the measured frequencies were less than the theoretical collsionless frequencies. This shift in frequency is then interpreted to be due to the fact that the system is shifted from the collisionless regime towards the hydrodynamic regime. Collisional damping associated with this shift was calculated and compared to the experimental damping. It is clear from Fig. 4 that below the transition temperature this kind of explanation is not possible, since most of the experimental points are not located between the collisionless and hydrodynamic curves.

To investigate the possibility that this damping might be due to the condensate-noncondensate atomic collisions, we now study the effect of these collisions by taking $`\alpha `$ to be nonzero. For nonzero $`\alpha `$ the dispersion relation takes the general form

$`\left({\displaystyle \underset{k=0}{\overset{9}{}}}P_C^{(k)}\alpha ^k\right)(\omega \tau _{22})^2`$ $`+`$ $`i\left({\displaystyle \underset{k=0}{\overset{8}{}}}P_I^{(k)}\alpha ^k\right)\omega \tau _{22}`$ (50)
$``$ $`{\displaystyle \underset{k=0}{\overset{7}{}}}P_H^{(k)}\alpha ^k=0.`$ (51)
We note that the $`k=0`$ term in this expression corresponds to Eq. (46), where $`\alpha `$ was indeed taken to be zero. Frequencies and damping rates can now be obtained by writing $`\omega `$ in terms of a real and imaginary parts, $`\omega =\omega _r+i\omega _i`$, and then inserting this expression in the last equation. In general, one obtains expressions for $`\omega _r`$ and $`\omega _i`$ in terms of $`\overline{\omega }\tau _{22}`$ and $`\alpha `$.

We start our discussion of the frequencies and damping rates for $`T/T_{BEC}=0.79`$, since at this temperature the experiment provides data for both the condensate and noncondensate oscillations. We present the results of our calculation for the $`m=0`$ and $`m=2`$ modes of the thermal cloud in Fig. 5, where we have eliminated $`\overline{\omega }\tau _{22}`$ and plotted $`\omega _i`$ as a function of $`\omega _r`$ for $`\alpha =0`$ and $`\alpha =0.25`$. On the same plot we show the two experimental points corresponding to the $`m=0`$ and $`m=2`$ modes of the noncondensate oscillation. In this figure, the right end of the curves represent the collisionless regime ($`\overline{\omega }\tau _{22}1`$) and the left end is the hydrodynamic regime ($`\overline{\omega }\tau _{22}1`$). Therefore, the location where the dashed curves intersect with the horizontal axis can also be seen in Fig. 4 at $`T=0.79T_{BEC}`$. For the experimental point corresponding to the monopole mode, the error bars do not intersect with the dashed curve ($`\alpha =0`$), but they do with the solid curve ($`\alpha =0.25`$) indicating that a nonzero value of $`\alpha `$ is needed to account for the experimental damping. This indicates that condensate-noncondensate atomic collisions are mainly responsible for the observed damping of this mode. This is less clear for the quadrupole mode although the experimental data is certainly not inconsistent with this conclusion.

Next we show in Fig. 6 at the same temperature both the condensate as well as the noncondensate $`m=0`$ and $`m=2`$ modes for values of $`\alpha `$ ranging from 0 to 0.25. The effect of $`\alpha `$ on the noncondensate modes is a slight shift and rotation of the curves but their general structure is roughly preserved. This can be seen more clearly in Fig. 5. The effect of $`\alpha `$ on the condensate modes is much larger, it gives rise to a large upward shift for the whole curve. The effect of $`\overline{\omega }\tau _{22}`$ on these curves becomes smaller for larger values of $`\alpha `$, which can be seen by noting that the radius of the small semicircles becomes smaller for larger values of $`\alpha `$. Again, we see that a nonzero value of $`\alpha `$ provides sufficient damping to account for the experimental observations. However, there is a discrepancy between our calculation and the experiment in the value of the frequency of the $`m=0`$ mode, since our calculation predicts $`\omega /\omega _11.68\pm 0.4`$, whereas the experimental data give a frequency $`2.05\omega _1`$. This discrepancy was also present in the calculations of Refs.  and is discussed in Ref. . The authors of the last reference suggested that a possible reason for this discrepancy may be that the observed mode is in fact an in-phase $`m=0`$ mode (the upper solid curve in Fig. 4) rather than an out-of-phase one. This suggestion was based on calculating the oscillator strenghts of the in-phase and out-of-phase $`m=0`$ modes. It turned out that it is indeed possible experimentally to excite both modes simultaneously in the temperature range $`T0.25T_{BEC}`$ to $`T0.5T_{BEC}`$. Although the calculation in Ref.  was performed in the collisionless limit, where there is no damping, we expect that with damping this argument remains qualitatively correct.

For the rest of the experimental points we can, in principle, explain the data using a two-parameter fit, namely with $`\alpha `$ and $`\overline{\omega }\tau _{22}`$. In fact, the collisions between atoms from the condensate with those from the thermal cloud are the main cause of damping in the condensate modes, and we can even explain the data using a one-parameter fit, namely only $`\alpha `$, by assuming that the system is in the collisionless regime with respect to $`\overline{\omega }\tau _{22}`$, i.e., $`\overline{\omega }\tau _{22}1`$. We perform this calculation by using a function $`\alpha (T/T_{BEC})`$ that gives the best fit for the experimental damping rates. In Fig. 7 we plot the function $`\alpha (T/T_{BEC})`$ and the resulting temperature-dependent frequencies and damping rates. This figure contains the main results of this paper. First of all we notice that now, with a nonzero $`\alpha `$, we obtain better agreement with the experimental data for the quadrupole mode frequencies than before where $`\alpha `$ was taken to be zero. This can be seen by comparing this figure with Fig. 4. We note also that, depending on the value of $`\alpha `$, the mode frequency may shift upwards or downwards. By comparing Fig. 4 with Fig. 7 for the mode frequencies one can clearly see that for temperatures below approximatelly $`0.6T_{BEC}`$ the theoretical curves are shifted upwards, whereas for the higher temperatures the curves are shifted downwards, in the end giving rise to the good agreement of the quadrupole mode with experiment. Secondly, we notice that the same function $`\alpha (T/T_{BEC})`$ also gives a good fit for the damping rates of both the monopole and quadrupole modes. Finally, there is still a discrepancy between our predictions and the experimental findings for the frequencies of the monopole mode at higher temperatures, but this may be resolved by a reinterpretation of the observed mode, as explained previously.

Finally, we want to mention an interesting feature of the in-phase and out-of-phase $`m=0`$ modes which can already be seen in Fig. 4 at $`T/T_{BEC}0.05`$. It is the familiar anticrossing tendency of these modes obtained also by Bijlsma and Stoof and Zaremba et al. . It is more dramatic here than in these two works due to the presence of the new parameter $`\alpha `$. This behavior can be seen in Fig. 8 where we see the two modes come close to each other and eventually intersect at some ‘critical’ values of $`\alpha `$ twice.

## IV Conclusion

We have extended the method of calculating collisional damping, used previously above the Bose-Einstein condensation transition temperature , to below the transition temperature. Furthermore, we have included the effect of noncondensate-condensate collisions. By comparing with experiment, we conclude that it is presumably this collision process which is mainly responsible for the observed damping. Our theory provides a general dispersion relation that gives complex mode frequencies at a certain temperature as a function of two dimensionless parameters, namely $`1/\overline{\omega }\tau _{22}`$ and $`\alpha `$, that characterize the noncondensate-noncondensate and noncondensate-condensate atomic collisions, respectively. Our results for the fully collisionless frequencies agree with those of Ref.  for most of the temperature range below the transition temperature. The $`m=0`$ and $`m=2`$ hydrodynamic frequencies are, to the best of our knowledge, calculated here below the transition temperature for the first time.

At present we have not carried out a microscopic calculation of the mean collision times $`\tau _{22}`$ and $`\tau _{12}\alpha ^1`$, which are treated as phenomenological parameters here with the possibility to investigate the intermediate regime with respect to either of the two collision processes. However, the homogeneous calculations of Zaremba, Nikuni, and Griffin  indicate that the function $`\alpha (T/T_{BEC})`$ that we have obtained from a fit to the experimental data, has the correct order of magnitude and qualitatively also the correct temperature-dependence. Nevertheless such a microscopic calculation for the inhomogeneous experimental conditions of interest would be very desirable and is left for future work.

We obtain a rather good agreement with the experimental results for the mode frequencies and damping rates apart from the discrepancy for the $`m=0`$ out-of-phase mode for $`T>0.7T_{BEC}`$. While the experiment shows an upward shift in the frequency with respect to its zero temperature value, we predict a downward shift. Resolution of this discrepancy seems to require more accurate experimental data, which will also provide decisive comparison with the theory presented here.

## Acknowledgements

The authors would like to thank Michiel Bijlsma, Chris Pethick and Henrik Smith for useful suggestions and remarks.