Landau damping in trapped Bose condensed gases

We study Landau damping in dilute Bose–Einstein condensed gases in both spherical and prolate ellipsoidal harmonic traps. We solve the Bogoliubov equations for the mode spectrum in both of these cases, and calculate the damping by summing over transitions between excited quasiparticle states. The results for the spherical case are compared to those obtained in the Hartree–Fock (HF) approximation, where the excitations take on a single-particle character, and excellent agreement between the two approaches is found. We have also taken the semiclassical limit of the HF approximation and obtain a novel expression for the Landau damping rate involving the time-dependent self-diffusion function of the thermal cloud. As a final approach, we study the decay of a condensate mode by making use of dynamical simulations in which both the condensate and thermal cloud are evolved explicitly as a function of time. A detailed comparison of all these methods over a wide range of sample sizes and trap geometries is presented.


Introduction
One of the more challenging problems in the study of Bose-Einstein condensation (BEC) in trapped atomic gases concerns the finite-temperature dynamics of a condensate interacting with a noncondensed, thermal component. Of particular interest in this regard is the study of collective oscillations of the condensate, where one would expect to see the influence of the thermal cloud on both the frequency and damping of the modes. Experiments which address these properties have been performed in a number of laboratories around the world [1]- [4], and provide an ideal test-bed for the development of theories of Bose condensed gases at finite temperatures.
One theoretical approach that can be used in the analysis of these experiments is the ZNG theory [5]. In this approach, the condensate is described by means of a generalized Gross-Pitaevskii equation while the dynamics of the thermal cloud is determined by a Boltzmann-like kinetic equation. The full numerical implementation of this theory has recently been presented [6] and its use in the analysis of a variety of experiments has appeared in a series of papers [7]- [10]. An important part of this previous work involved a detailed study of the collisional dynamics which is responsible for the equilibration of the system with respect to both energy and particle number. In this paper, however, we shall focus on an aspect for which collisions are of secondary importance, namely Landau damping. This damping mechanism is dominant in the low-density systems studied experimentally [1,3,4]. In this so-called collisionless regime, meanfield interactions mediate the transfer of energy from the condensate to the thermal component, leading to the damping of condensate collective modes.
The subject of Landau damping in dilute BECs has been explored by several authors. Early work on a uniform gas was conducted by Hohenberg and Martin [11] for low temperatures, while results for higher temperatures were obtained by Szépfalusy and Kondor [12]. Liu [13] used results obtained for the uniform gas to estimate the damping in the case of a trapped gas, while Pitaevskii and Stringari [14] and Fedichev et al [15] developed expressions for the Landau damping rate that explicitly accounted for the trap geometry. Similar results were derived by 88. 3 Giorgini [16,17] who also considered the Baliaev damping process whereby a collective mode decays into two lower-energy excitations. The latter can be ignored, however, for the low-lying collective modes in a trapped gas because of energy conservation restrictions. Semiclassical (sc) results for Landau damping were also obtained using a dielectric formalism by Reidl et al [18], who found quite good agreement with the results of experiment [1]. The ZNG theory also provides a description of Landau damping and applications of the theory to a number of different condensate collective modes [7]- [9] found damping rates that agreed well with experiment.
The perturbation theory approach of [14,16,17] provides an expression for the Landau damping rate that involves a sum over transitions between pairs of Bogoliubov excitations. The first fully quantum mechanical evaluation of this expression was performed by Guilleumas and Pitaevskii (which we shall henceforth denote as GP I) [19], for a spherically symmetric trapping potential. We should also mention that a related calculation was performed by Das and Bergeman [20], obtaining similar results to GP I. In the ZNG theory, the description of Landau damping is quite different since it involves the dynamical evolution of the condensate in the presence of a thermal cloud treated semiclassically at the level of the Hartree-Fock (HF) approximation. The relation of this approach to that used in GP I is far from evident, although it is clear that the physical bases of the two are the same. To investigate this question, we performed numerical simulations in an earlier paper [6] for the system studied in GP I and in fact found quantitative agreement. One of the aims of the present paper is to compare the two methods over a much larger range of condensate sizes, which means that we have had to repeat calculations of the kind performed in GP I. By doing so, we are then able to investigate the consequences of using HF, as opposed to Bogoliubov, excitations in the Landau damping calculation. A reformulation in terms of HF excitations has the added benefit of allowing the sc limit to be taken in a straightforward way, thereby completing the connection to the ZNG theory.
Perhaps an even more interesting question concerns the damping rate in anisotropic traps. This issue arose in experiments on the transverse breathing mode in a highly elongated trapped gas [4], where the damping rate was found to be anomalously low-around an order of magnitude smaller than one would expect in similar circumstances on the basis of both experiment and theory. In [9] we simulated this experiment and found that the low damping rate could be explained as a consequence of an in-phase collective oscillation of both the condensate and thermal cloud. On the other hand, if the thermal cloud is not set into oscillation, the Landau damping rate is found to have a value typical of those seen in other experiments. This result apparently contradicts the conclusion of Guilleumas and Pitaevskii (GP II) [21] who performed a Landau damping calculation for an infinitely long cylindrical condensate. They found a damping rate which was extremely low, around an order of magnitude smaller than the experimental results, and two orders of magnitude smaller than our simulation results for a stationary thermal cloud. A natural question concerns the validity of modelling the experimental geometry by an infinite cylindrical trap which changes the nature of the excitation spectrum and presumably the damping rate. In this paper we address this question by repeating the GP I calculations as a function of the trap anisotropy and compare these results with those obtained from our simulations. Our tentative conclusion is that the damping rate has an anomalous dependence on anisotropy in the limit of the infinite cylinder, so while the GP II result is correct, it does not apply at the experimental anisotropy. 88.4

Landau damping of Bogoliubov excitations
For purposes of completeness it is useful to summarize the Bogoliubov theory and the way in which a condensate collective mode is damped by thermal excitations. In a second-quantized notation, the Hamiltonian is given bŷ where V ext (r) is the trapping potential and g = 4π ah 2 /m is the interaction parameter defined in terms of the s-wave scattering length a. If the field operator is expressed aŝ where 0 is the condensate wavefunction andψ is the excitation (or fluctuation) field operator, the Bogoliubov approximation is generated by expandingĤ to second order inψ. The terms linear inψ vanish if 0 satisfies the Gross-Pitaevskii equation where n c0 = | 0 | 2 is the condensate density. At this level of approximation, the condensate does not interact directly with the thermal excitations as determined in the Bogoliubov theory. The resulting Bogoliubov Hamiltonian is then diagonalized by means of the Bogoliubov transformationψ The quasiparticle amplitudes u i and v i satisfy the Bogoliubov equationŝ where we introduce the operatorL ≡ −h 2 ∇ 2 /2m + V ext + 2gn c0 − µ. The orthonormality of the quasiparticle amplitudes is specified by the relation Apart from a constant, the Bogoliubov Hamiltonian is given bŷ where the excitation energies E i correspond to condensate collective modes. In thermal equilibrium, these modes are populated according to a Bose distribution The cubic terms in the expansion ofĤ , couple the Bogoliubov excitations and lead to a mechanism for their decay. In particular, a mode which has been excited, and therefore is no longer in equilibrium with the other thermal 88.5 excitations, will experience a decay known as Landau damping. To represent this situation, one defines a state in which the mode occupation n osc is large compared to the equilibrium value. The energy in this mode, E osc =hω osc n osc , then decays as n osc goes to equilibrium. The rate of transition from the initial nonequilibrium state to any other state can be determined by perturbation theory [14] and the average rate of change of the energy in this mode is found to bė where the transition matrix element is For future reference we note that the HF approximation is recovered by setting v i (but not v osc ) in (10) equal to zero and determining u i from the first of the Bogoliubov equations in (5) with The damping rate is defined according to 2 = −Ė osc /E osc (which implies that the amplitude of the mode decays as e − t ). Thus, This expression is the basis of the calculation of Landau damping as carried out by Guilleumas and Pitaevskii for both spherically symmetric and cylindrical traps [19,21]. In order to display the results of our calculations it is convenient to follow their notation: where is a 'damping strength' associated with each possible transition, with frequency difference The calculation of the damping then consists of solving the Bogoliubov equations (5) for u osc , v osc and ω osc for the classical mode of interest, together with the corresponding quantities for the thermally populated excitations. The matrix elements A i j can then be evaluated for transitions between each pair of excitations, which in turn yields γ i j . More details of these calculations will be given in the next section, while the calculation of the total damping rate , including how we deal with the delta function in (12), will be presented in section 3.

Anisotropic traps
In performing the above calculations, a significant saving in computational effort can be realized by considering a spherically symmetric trap. Not only does (10) effectively reduce to a onedimensional integral, but the spherical symmetry also reduces the number of excitations that have to be considered. In particular, if we label the modes with the usual quantum numbers {n, l, m}, then states corresponding to a given l are (2l + 1)-fold degenerate, and A i j need only be calculated explicitly for one of these states. This spherical case was first studied in [19], and we will repeat some of these calculations in section 3.1 in order to compare various 88.6 approximations. However, here we are also interested in the damping in anisotropic traps which is the most common experimental situation. It turns out that calculations for anisotropic traps, even those with axial symmetry, are much more involved, since the degeneracy of excitations with different m quantum numbers is lifted. As a result, determining the Bogoliubov excitation spectrum requires the diagonalization of much larger matrices than in the spherical case. The corresponding increase in the number of possible transitions, along with the fact that evaluation of the matrix elements A i j now involves integrating over both radial and angular variables, significantly increases computational overheads. In the remaining part of this section we shall describe how we tackle these problems by making use of a spherical harmonic basis. The calculation for an isotropic trap is then a simplified version of this general scheme.
As a first step, we find the ground state condensate wavefunction as detailed in [22]. This makes use of a set of normalized basis functions ψ nlm (r) = R nl (r )Y lm (θ, φ), which are eigenfunctions of the Hamiltonianĥ 0 = −h 2 ∇ 2 /2m + V 0 (r ). The potential V 0 (r ) is the l = 0 component of the total axially symmetric effective potential V (r) ≡ V ext (r) + gn c (r) when expanded in terms of Legendre polynomials, V (r, θ) = l V l (r )P l (cos θ). We take the external harmonic potential to be of the form V ext (r) = (mω 2 ⊥ /2)(x 2 + y 2 +λ 2 z 2 ), where λ is the anisotropy parameter. By defining a nonspherical perturbation V (r) = V (r) − V 0 (r ), and expanding an arbitrary solution of the GP equation as φ(r) = nl a nl ψ nlm (r), one can set up a matrix equation for the coefficients a nl . The condensate wavefunction is then given by the lowest-energy evenparity solution 0 (r) = √ N c φ 0 (r), which in turn can be used to calculate the condensate density n c (r). Updating the effective potential V (r) then leads to a new matrix equation, which is solved to yield an updated wavefunction. The calculation is repeated until 0 converges to the correct condensate wavefunction, with the lowest eigenvalue giving the chemical potential µ. In the process, one also generates a full set of eigenfunctions φ α (r) for the Hamiltonian Once we have this information, the Bogoliubov equations can be solved by introducing . The index i labels the different excitations and includes the azimuthal quantum number m. Again, this leads to a matrix equation, which can be solved to obtain the eigenfunctions ψ ± i (r) and energy eigenvalue E i for each mode. For the purposes of the Landau damping calculation, which we move onto next, it is convenient to express these functions as by making use of the expansions of φ α in terms of the ψ nlm . For the particular mode of interest (i.e., the one for which we are calculating the damping) we can define the mode densities, δn + (r) = 0 (r)ψ + osc (r) and δn − (r) = 0 (r)ψ − osc (r). Using this notation, the Landau damping matrix elements become In evaluating this integral, it is convenient to expand the mode densities (where the mode has a particular m-number, which we denote asm) in terms of spherical harmonics: Substituting (14), (15) and (17) into (16) gives where in which l 1 l 2 m 1 m 2 |l 3 m 3 is the Clebsch-Gordan coefficient [23], and where m j = m i +m. The calculation of the Landau damping then involves evaluation of (19) for each pair of excitations with a difference of energies in a particular interval E min < E j − E i < E max . In this paper we shall focus on the Landau damping of modes withm = 0 and even parity, so selection rules connect excitations with m = 0 and the same parity. Unlike in the spherically symmetric case, however, excitations with different m-values must be summed explicitly. Nevertheless, since the contributions from the ±m i excitations are identical, a twofold saving can be gained for m i = 0.
Our results for the breathing mode as a function of anisotropy will be presented in section 3.2.

Hartree-Fock approximation
In this section we obtain an expression for Landau damping in the HF approximation. Although this can be done trivially by setting v i = 0 in (10), it is useful to rederive the result using an alternative method for a number of reasons. First it clearly shows that Landau damping is associated with the work done on the thermal cloud by the dynamic mean field of the oscillating condensate. Second, this reformulation allows one to take the sc limit in a straightforward way, thereby facilitating a comparison with the sc simulations that we have performed on the basis of the ZNG theory. Within the HF approximation, the condensate and thermal cloud are treated as two distinct components. The condensate evolves dynamically according to a time-dependent GP equation while the thermal cloud responds to the condensate mean field. This picture is quite distinct from the Bogoliubov approach in which the 'condensate mode' is distinguished from the 'thermal excitations' only through the assumed different occupations of the respective modes. Common to both pictures, however, is the absence of a thermal cloud mean field acting back on the condensate. Such effects are included in the ZNG theory and are crucial for the satisfaction of the generalized Kohn theorem [5].
For small deviations from equilibrium, the time-dependent condensate wavefunction is given by where the fluctuation δ is obtained from the linearized GP equation. The fluctuating part of the condensate density is then 88.8 which gives rise to a dynamic mean field 2g δn c (r, t) acting on the thermal cloud. Thus the thermal cloud experiences a perturbation wheren(r) represents the density operator of the thermal component. The linearized GP equation has solutions of the form where the u and v are, apart from a normalization, the same Bogoliubov amplitudes as introduced earlier and ω osc is the frequency of the particular condensate mode of interest. In terms of this wavefunction, the condensate density fluctuation is given by With this form of the condensate density fluctuation, the perturbation given in (23) has a harmonic time dependence and can be treated by means of time-dependent perturbation theory. One thus finds that the time-averaged rate of change of the thermal cloud energy is given to lowest order in the perturbation by the expression where χ (r, r , ω) is the imaginary part of the time Fourier transform of the thermal cloud density response function The angular brackets, · · · 0 , denote an equilibrium expectation value. Treating the thermal atoms as independent HF excitations, the spectral density is given by the expression where φ i (r) is the HF wavefunction (as obtained from (5) by setting v i = 0) for the i th state with energy ε i , and f i is the thermal Bose occupation number. Substituting this expression into (27), with As stated earlier, this matrix element can be obtained from (10) by simply setting v i = 0, but retaining v osc which is needed to properly define the condensate density fluctuation. The rate of change of the condensate energyĖ osc is of course the negative of (30). To extract a damping rate we must divide this by the energy of the mode which is determined by the amplitude of the condensate density fluctuation. If u and v are normalized according to (6), this energy is justhω osc . 88.9

Semiclassical HF
Our purpose in this section is to show how the quantal HF damping can be reduced to a sc form. This will allow us to make contact with the sc ZNG simulations. The response function in (28) can be written as This correlation function can be expressed in the form [24] φ(r, which allows the semiclassical limit for the dynamics to be taken straightforwardly. Settinḡ where now n(r, t) is a classical density variable. The Fourier transform of (34) gives To evaluate the correlation function in (35), we consider a thermal cloud containingÑ atoms. The density of the cloud at time t is thus given by where r i (t) is the position of the i th particle at time t. Thus, Since the atoms in the cloud are independent and equivalent, we have where we have defined the self-diffusion function 2 [25] andñ 0 (r) is the equilibrium thermal cloud density. Since the product of equilibrium densities is time independent, it will not contribute to the Fourier transform in (35) at finite ω and can be dropped. Combining these results, the sc Landau damping rate is given by wheref (ω) is the Fourier transform of the function 88. 10 We thus see that the quantity determining the damping rate is an appropriately weighted spatial average of the self-diffusion function G s . The latter can be evaluated analytically for some simple model situations which provide insight into its behaviour. However, in the present situation, the thermal atoms are moving in the combined potential of the trap and the mean field of the condensate and G s can only be determined numerically. It is then easier to deal with f (t) directly.
Using the definition of G s , f (t) can be expressed as which is convenient for numerical evaluation. We note that r 1 (t) is the position of an atom at time t starting from some initial position r 1 (0) ≡ r 0 with momentum p 1 (0) ≡ p 0 . The initial phase point is distributed according to the equilibrium Bose distribution f 0 (r 0 , p 0 ), normalized toÑ , and so the correlation function takes the form Thus, f (t) can be determined by averaging over all possible trajectories starting from some initial phase point. In practice this can be done by considering an ensemble of N T test particles distributed in phase space according to f 0 (r 0 , p 0 ) and then performing the discrete summation N T is chosen sufficiently large to ensure that statistical fluctuations in f (t) are acceptably small. Examples of such calculations will be given later. Finally, a numerical Fourier transform of f (t) at ω = ω osc yields the semiclassical damping rate.
The value of f (0) can also be evaluated directly as This serves as a check of the sum over test particles described above. In addition, we note that f (t) tends to a finite limiting value for t → ∞. This limiting value can be subtracted from f (t) in the evaluation off (ω) since this only affectsf (ω) at ω = 0.

ZNG simulations
Our final method for evaluating the Landau damping is based on N -body simulations [6] in which collisions are neglected. In this numerical scheme, the condensate wavefunction is evolved in time by numerically solving the generalized GP equation using a split-operator FFT method, while the evolution of the thermal cloud is calculated by sampling with a large number of test particles, which are evolved classically. This is in fact equivalent to solving a collisionless Boltzmann kinetic equation for the thermal cloud, where the thermal excitations are described according to the sc HF approximation (i.e., as single-particle excitations in a quasi-uniform gas [5]). The dynamic mean field of the thermal cloud then gives rise to decay of the condensate oscillation, which is identified with Landau damping. Full details of the numerical methods are given in [6]. The first step in the calculation is the self-consistent determination of the condensate wavefunction together with the thermal cloud density. The simulation is then initiated by exciting 88.11 the condensate in such a way as to realize the relevant mode of study. It is straightforward to show that this can be achieved by adding ψ − osc (r), as obtained from the solution of the Bogoliubov equations (5), to the condensate ground state wavefunction. The subsequent evolution of the mode of interest can then be monitored by evaluating the moments x 2 , y 2 and z 2 of the condensate as a function of time. These are calculated numerically using χ (t) = dr χ | (r, t)| 2 . A plot against time yields a damped oscillation, which is fitted to a sinusoidal function with an exponentially decaying factor exp(− t). Each simulation extends for a time of ω ⊥ t = 10, where the relatively short time is chosen so as to avoid driving the thermal cloud too far out of equilibrium. This point is discussed more fully in [6]. The damping rate, , is then extracted from the simulation data, which can then be compared to the results of our other calculations. It should be emphasized that within this approach the damping of the condensate oscillation is a direct result of the dynamic mean field that the thermal cloud exerts on the condensate. In other words, no damping would occur if the condensate were to simply oscillate in the static equilibrium mean field of the thermal cloud.

Isotropic traps
We now move on to describing the results of our calculations, where we begin with the case of the Bogoliubov method for an isotropic trap (λ = 1, ω 0 ≡ ω ⊥ = 2π × 187 Hz). The specific system being considered is a gas of 87 Rb atoms with a scattering length a = 5.82 × 10 −9 m. The calculation follows closely that of Guilleumas and Pitaevskii (GP) [19]. The results are conveniently displayed as a histogram which plots the damping strengths, γ i j , versus the excitation frequency ω i j = (E j − E i )/h. Figure 1(a) shows such a plot for N c = 2.5 × 10 5 condensate atoms at a temperature of k B T /µ = 1.5, where the chemical potential takes its Thomas-Fermi value µ = 29.9hω 0 . The spectrum features several high peaks, which correspond to transitions between low-lying excitations with large overlaps (i.e., large matrix elements A i j ) together with large population factors, f i − f j . Following GP, we shall refer to these peaks as 'resonances'. They are not expected to play a significant role in the Landau damping process, unless one happens to lie very close to the mode frequency ω osc . For this reason we shall ignore these transitions when calculating the Landau damping, and instead focus on the small-amplitude quasi-continuous 'background'.
A conceptual difficulty arises in the calculation in that the excitation spectrum in (11) consists of discrete delta functions. Taken literally, this would imply that the Landau damping rate would either be zero or infinite depending on the location of the oscillation frequency ω osc . Following [19], this difficulty is overcome by replacing each delta function in the background spectrum by a Lorentzian /{(2π)[(ω i j − ω osc ) 2 + 2 /4]}. Physically, this accounts for the fact that the excitations will themselves have a finite lifetime in a more rigorous theory. As a result, the Landau damping rate is determined by averaging over many peaks weighted according to their proximity to the mode frequency. The width factor, , however is somewhat arbitrary, and as shown in figure 2, the result for will vary with . Nevertheless, one can see that the variation is weak when /ω 0 lies between 0.05 and 0.20. We fit the data in this range to a straight line, and extrapolate back to = 0 to yield an estimate of the Landau damping rate. The result as a function of temperature is shown in figure 3, for different numbers of condensate atoms. One sees in the plots the typical rapid increase of the rate at low temperatures followed by a nearly linear temperature dependence at the higher temperatures. Our results for N c = 5 × 10 4 are in complete agreement with those of Guilleumas and Pitaevskii [19].
We have repeated the calculation for HF thermal excitations following the method discussed in section 2.3. The spectral density in this case is shown in figure 1(b), for the same parameters as in (a). The main difference between the two plots is the absence of the strong resonances seen in the Bogoliubov calculation. However, the background spectrum is remarkably similar, especially in the vicinity of the mode frequency ω osc . Insight into why this is happening is provided by comparing the Bogoliubov and HF excitation frequencies, as plotted in figure 4 for N c = 5 × 10 4 . As found previously [26], the two spectra are completely different at low energies and angular momenta, but converge at high E and l, demonstrating that the excitations take on a single-particle character [27,28] in this region. It is these excitations that are responsible for the majority of the background transitions in figure 1 which explains the similarity between the HF (calculated in the same way as described above) and Bogoliubov damping rates plotted in figure 3. The excellent agreement with the Bogoliubov results for a wide range of condensate sizes is in contrast with the case for the uniform Bose gas [14], where one finds a significant difference between the two approaches. The reason for this is that the Bogoliubov spectrum for a uniform gas only approaches the single-particle HF form at temperatures much larger than the chemical potential, while in a trapped gas, surface excitations with high multipolarities are also important, even at low temperatures. We thus arrive at the important conclusion that the 'collective nature' of the excitations has little bearing on the determination of Landau damping in trapped Bose gases.
It is also of interest to compare the Bogoliubov results to those of our sc simulations based on the ZNG theory (as described in section 2.5). We plot the damping rates for both calculations  as a function of the number of condensate atoms in figure 5, for three different temperatures. The behaviours seen for the two sets of calculations are similar, with a more rapid rate of increase of the damping rate at low N c followed by a less rapid increase at higher N c . The data at N c = 5×10 4 and 1.5×10 5 correspond to those plotted in figure 8 of [6], where we concluded that the two approaches were in good agreement. The comparison in figure 5 shows that the agreement persists over a wider range of condensate number, although some systematic differences do exist. The simulation rates at low N c tend to be lower than the Bogoliubov damping rate, while at high N c there is a tendency for them to be slightly larger. These small differences are presumably due to the different approximations used in the two sets of calculations. These include (i) the use of HF as opposed to Bogoliubov excitations in the simulations, (ii) the use of the sc approximation for the thermal cloud dynamics and (iii) the inclusion of the thermal cloud mean field in the selfconsistent calculation of both the equilibrium and dynamical properties. Concerning (i), we saw in figure 3 that there is little difference between the Bogoliubov and HF results over a wide range of temperatures when the excitations are both treated quantum mechanically. We therefore do not expect the use of Bogoliubov, as opposed to HF, quasiparticles to appreciably change the results.
The effect of the sc approximation itself, however, is less clear. To isolate this effect we make use of the formulation outlined in section 2.4. The calculation is based on equation (40) with f (t) determined from (44) by sampling the mode density δn − (r) along classical trajectories. To be specific, we again consider the radial breathing mode in an isotropic trap with frequency ω 0 = 2π × 187 s −1 . The mode density being sampled in this case is shown in figure 6 for N c = 10 5 . The node in the mode density is of course required in order for the volume integral of the mode density to vanish.
We consider first a model in which the thermal atoms move in a purely harmonic potential; that is, we ignore the mean field of the condensate. Furthermore, we assume that the phase space density of the atoms is given by a Boltzmann distribution. With these assumptions, the self-diffusion function G s (r, r , t) can be evaluated analytically with the result [29] This function is periodic with period T 0 = 2π/ω 0 : it recovers its t = 0 value G s (r, r , 0) = δ(r − r )ñ 0 (r) at the times t n = nT 0 . This is a special property of purely harmonic confinement. When substituted into (41), the resulting function f (t) has period T 0 /2 since the mode density δn − (r) has inversion symmetry. We show in figure 7(a) f (t) as calculated from (44) by following   of the dominant 2ω 0 frequency of the purely harmonic case. The mean field of the condensate disrupts the periodicity of the particle orbits that occurs for harmonic confinement and, as a result, f (t) saturates at long times to a constant limiting value. The inset of figure 7(b) shows f (t) on a smaller timescale; the behaviour seen is consistent with f (t) being an even function of time.
In figure 8(a) we show f (t) for the same conditions as in figure 7(b), but with the Boltzmann distribution replaced by the Bose distribution. The behaviour of f (t) is qualitatively very similar; the main difference is the limiting long-time value which can be attributed to the different thermal distributions in the two simulations. It is worth commenting on what this asymptotic value is due to. If the position of an atom at long times is uncorrelated with its starting point, one might expect the self-diffusion function to tend towardsñ 0 (r)ñ 0 (r )/Ñ . This would yield an 'equilibrium' value of f (t) of We find, however, that f (t) in our simulations does not tend to this limit. The reason for this is that long-time correlations persist since a given atom retains its initial energy in the course of the dynamical evolution. Even if the nonharmonic perturbation were to lead to an ergodic distribution on a given energy surface, the equilibrium limiting form can only arise if equilibrating collisions between thermal atoms take place. This was confirmed by performing a simulation in which collisions are included following the methods discussed in [6]. The dashed curve in figure 8(a) shows the result of this calculation and we now indeed find that f (t) → f eq as t → ∞. In figure 8(b) we show the Fourier transform of f (t). It has peaks at frequencies close to multiples of 2ω 0 as would be expected in view of figure 8(a). Less obvious is the fact thatf (ω) appears to be positive definite. This behaviour, however, is necessary sincef (ω) is just the sc limit of the spectral density in (30). It is therefore reassuring that this property emerges from our simulations. The dashed curve in figure 8(b) is the corresponding Fourier transform when collisions are included. The main difference between the two curves occurs at low frequencies where the collisional curve has a Lorentzian peak due to the relaxation of f (t) to f eq on a timescale of ω 0 t ∼ 10.

88.17
The sc Landau damping rate is determined by the value off (ω) at ω = ω osc for the radial breathing mode. This frequency is indicated by the vertical line in figure 8(b) and it is clear that collisions are not important in this case. The Landau damping rates calculated on the basis of (40) are plotted as a function of temperature in figure 9 which also displays the Bogoliubov results for comparison. The two sets of calculations are qualitatively similar and show a similar dependence on the condensate number N c . However, the sc results underestimate the Bogoliubov results at higher temperatures.
To understand these differences we show in figure 10 excitation spectral densities for the three approximations over a much larger frequency range. For the Bogoliubov and HF calculations, the spectral densities are given by (11), with ω osc replaced by an arbitrary frequency ω, and convolved with a Lorentzian of width = 0.1ω 0 . This representation is more informative than the histograms in figure 1 since it reveals the total spectral weight as a function of frequency.  The corresponding sc spectral density follows from (40) and is given by It is clear from this figure that the spectral densities in all three approximations are qualitatively similar, with a series of broad peaks near multiples of 2ω 0 . However, both the Bogoliubov and HF calculations have a higher spectral density in the range 2ω 0 < ω < 4ω 0 than the sc result, which accounts for the smaller sc damping rates shown in figure 9. The larger quantum spectral 88.19 densities are associated with low-energy states for which the zero-point energy and the effects of barrier penetration are more important. These quantum effects are of course missed in the sc HF calculation. However, the overall similarity of the curves in figure 10 indicates that the sc HF approximation is providing a good qualitative description of the thermal cloud excitations.
In order to gain more insight into these results, we have performed another simulation which incorporates the same physics as the sc HF approximation formulated in section 2.4. As in the full simulations in section 2.5, we evolve the condensate using the time-dependent GP equation but do not include the mean field of the thermal cloud acting on the condensate. The condensate thus oscillates with a fixed amplitude giving rise to a harmonic perturbation of the thermal cloud which itself starts off as an equilibrium distribution and evolves classically in time in the presence of the dynamic mean field of the condensate. These are the same ingredients as enter the HF perturbation theory calculation. We then monitor the thermal cloud energy is the equilibrium effective potential acting on the thermal cloud. As a result of the work done by the condensate on the thermal cloud, the latter experiences a secular increase in energy which we can identify with the time-averaged rate of energy transfer in (30). Some points obtained in this way are shown in figure 9. These points agree quite well with the Bogoliubov rates, and by the same token, with the full simulations as shown in figure 5. Thus, the two types of simulation are mutually consistent but deviate from the results based on the sc HF expression in (40). Apparently, following the dynamics of the thermal cloud in the presence of an oscillating condensate is not entirely equivalent to the evaluation of the equilibrium correlation function in the sc HF expression. At the present time we do not have a satisfactory explanation for the observed difference. We speculate that it may be related to the nonequilibrium nature of the thermal cloud distribution in the simulations, but we have not been able to find a way to check this.

Anisotropic traps
We now describe the effect of introducing an anisotropy (λ < 1) such that the condensate takes on a prolate geometry. The anisotropy lifts the degeneracy with respect to the quantum number m which necessitates the inclusion of a much larger number of transitions in our Bogoliubov calculations. As a result, calculations for N c > 10 4 are numerically prohibitive. The same is true for large anisotropies, so we are limited here to 0.3 < λ ≤ 1. This range, however, is sufficiently large to reveal some general trends.
The effect of imposing a slight anisotropy on the excitation spectrum is shown in figure 11 (for N c = 10 4 and with the geometric mean of the trap frequency ω ho = ω ⊥ λ 1/3 kept fixed at 2π ×187 Hz). One sees that some of the larger resonances, corresponding to transitions between the lower-energy excitations, are each split into a series of shorter peaks. This is most clearly illustrated for the large peak at ω i j /ω ho = 1.94, which is associate with a pair of low-lying l = 1 modes. The anisotropy splits this peak into three since the m = −1, 0 and 1 modes are no longer degenerate. As the anisotropy increases, the peaks continue to separate, as shown by the bottom panel in figure 11 for λ = 0.9.
An obvious consequence of this is that the spectrum is now much denser, but one can extract the damping rate in much the same way as described above for isotropic traps. The only difference is that we no longer ignore the larger peaks (resonances) in our calculation, since these are more numerous and less distinguishable from the background. There is therefore no simple criterion that can be used to remove them. The damping rate for each λ is plotted in   figure 12, where one sees a downward trend as the anisotropy increases towards a cigar-shaped geometry (decreasing λ). The anomalously high damping rate at λ = 0.55 is a consequence of a mixing of the radial breathing mode with another m = 0 mode which happens to possess a similar frequency at this anisotropy. We show the effect of this mode crossing on the spectrum by scanning through λ = 0.55 in figure 13. Many more high resonances are seen at λ = 0.55 than to either side, which in turn leads to the much higher damping rate shown in figure 12. We have here a case in which the 'resonant modes' do in fact contribute to the Landau damping. The shift of the spectrum towards the right is due to the ratio ω ⊥ /ω ho increasing as λ increases.
We also believe that the scatter in the Bogoliubov data seen in figure 13 can be accounted for in terms of resonant peaks approaching or receding from the vicinity of the breathing mode frequency.
The straight line in figure 12 is a least-squares fit to the Bogoliubov data (excluding the high-lying point at λ = 0.55). Extrapolating the straight line fit to λ = 0 suggests a damping rate of λ/ω osc 0.006 in this limit, which is around a factor of 4 smaller than in the spherical case. For comparison, we have also plotted in figure 12 the results of sc ZNG simulations. These calculations are again carried out for the same condensate number (N c = 10 4 ) and temperature (k B T /µ = 1.5). We note from figure 5 that the simulation damping rate at λ = 1 lies below the Bogoliubov calculation for this low number of condensate atoms. This seems to impart a slightly different λ-dependence as compared to the Bogoliubov results, but we still see a downward trend for λ < 0.5. In particular, there is a rapid decrease in the simulation damping rate for λ → 0.
These results shed some light on differences found in previously reported calculations [9,21]. Guilleumas and Pitaevskii [21] performed a Bogoliubov calculation for a cylindrical condensate (corresponding to λ = 0) and found a Landau damping about two orders of magnitude smaller than the typical values in figure 5, and an order of magnitude smaller than the damping rates observed experimentally for λ = 0.065 [4]. On the other hand, we found [9] a conventional Landau damping rate from ZNG simulations (as discussed in this paper) that was consistent with figure 5, but an order of magnitude larger than that from experiment. The rapid variation of the simulation results in figure 12 for small λ is a possible resolution of these different theoretical predictions. If the damping rate continues to drop rapidly for λ → 0, one might find a result consistent with that of Guilleumas and Pitaevskii [21]. Unfortunately, we cannot push our calculations to lower λ in order to confirm this. It is nevertheless clear from all of our calculations that the Landau damping in the cylindrical geometry appears to be an anomalous limit.
As a final point, we should emphasize that the conventional Landau damping rate discussed here fails to account for the recent observations of the damping of the transverse breathing mode 88.22 in elongated condensates [4]. This could only be achieved by including the full dynamics of the thermal cloud in the ZNG simulation [9]. Conventional Landau damping deals with the oscillation of the condensate in an otherwise equilibrium thermal cloud; the condensate does work on the thermal cloud and loses energy to it. However, if the thermal cloud is itself set into motion as is the case in the experiments [4], a change in the damping rate is to be expected. The experimental observations [4] can only be accounted for if this effect is taken into account [9].

Conclusions
In this paper we have studied Landau damping of Bose-Einstein condensates in isotropic and anisotropic traps and have compared a variety of methods for calculating the damping. Calculations based on perturbation theory, involving sums over transitions between excited states, show excellent agreement between a Bogoliubov treatment of the excitations and the HF approximation over a wide range of temperatures and numbers of atoms. These results demonstrate that Landau damping in trapped Bose condensed gases is essentially a 'singleparticle' phenomenon in that the collective nature of the excitations is unimportant. A sc limit of the HF approximation has also been developed, leading to a novel expression for the damping rate in terms of a classical correlation function. A detailed comparison with the quantal HF results has shown that the excitation spectra in the two cases are qualitatively similar but differ at a quantitative level. In particular, the sc results tend to underestimate the quantum results at higher temperatures. We also explored damping by means of N -body simulations (as discussed in [6]) in two complementary ways. In the first, the condensate wavefunction is evolved using the time-dependent GP equation and its interaction with the thermal cloud leads to a decay of the condensate oscillation amplitude. In the second, the condensate oscillates with a fixed amplitude and the rate of increase of the thermal cloud energy is determined. The two methods give very similar results which agree with the quantum perturbation results over a wide range of temperatures and condensate sizes, reinforcing the conclusions of [6]. However, all of these results differ somewhat (see figure 9) from the sc HF formulation. At present we have no explanation for these differences.
We have also studied in detail the damping in cigar-shaped traps using the Bogoliubov approximation and sc simulations. These results show a decreasing trend in the damping as the trap becomes increasingly anisotropic. Over the range of λ in both sets of calculations, the damping decreases by about a factor of two, with the sc simulations indicating a more rapid decrease for λ → 0. This behaviour may reconcile the differences between the results obtained for highly elongated condensates [9] and those obtained [21] in the cylindrical limit (λ = 0).