Hydrodynamic collective modes for cold trapped gases

For a three-dimensional Bose gas the condition $na^{3}\ll 1$, where $n$ is the density, corresponds to a dilute and weakly interacting system. In this case, the mean-field result for the energy density at zero temperature is given by Bog1 ; PiSt1 ; PeSm1 ($\hbar=k_{B}=m=1$, see App. A for units)

Taking a derivative with respect to $n$, we get the chemical potential $\mu=g_{3D}n$ which enters the Gibbs-Duhem relation $\mbox{d}P=n\mbox{d}\mu$ for the pressure. The equation of state written in the grand canonical variables then becomes

which is of polytropic type $P\propto\mu^{\alpha+1}$ with $\alpha=1$. At zero temperature the system is completely superfluid and will be described by superfluid hydrodynamics.

The oscillation frequencies of this system in a spherical parabolic trap $V_{ext}=\frac{m}{2}\omega_{0}^{2}r^{2}$ are found by solving the eigenvalue problem

for the differential operator

acting on a function $g(z)$. Here $d=3$, see Sec. III. The operator $A$ depends on the equation of state through $P^{\mu}(\mu)$ and $P^{\mu\mu}(\mu)$, where a superscript denotes differentiation with respect to $\mu$. The former of these quantities corresponds to the density $P^{\mu}=n$ while their ratio $P^{\mu}/P^{\mu\mu}=\partial n/\partial P=c^{2}$ is related to the velocity of sound $c$. Since only this ratio appears in the operator $A$, the prefactor of Eq. (3) is of no importance and the frequencies will not depend on the interaction strength on the mean field level. Indeed, Stringari Str1 found the analytic expression $\omega_{nl}=(2n^{2}+2nl+3n+l)^{1/2}\omega_{0}$, where $n$ and $l$ are integers. Obviously, no thermodynamic conclusions can be drawn from this formula and it may at best help to determine the trapping frequency when interactions and thermal effects are known to be very small. Stringari’s formula is a special case of

which holds for $d$-dimensional spherically symmetric, harmonic traps and arbitrary polytropic index $\alpha$. As we will show in App. D, this formula can be obtained by applying a simple polynomial ansatz $g(z)=\sum_{k=0}^{n}a_{k}\bar{z}^{k}$ with $\bar{z}=z/R^{2}$ to Eq. (4) for $P\propto\mu^{\alpha+1}$. It is already known in the literature Fuc1 ; Hei1 .

Eq. (6) reveals that the possible oscillation frequencies of a trapped system form a discrete set. We will see in the following that this property is true for an arbitrary equation of state $P(\mu,T)$. This feature does not arise as a result of imposed boundary conditions but is related to the domain of definition of the operator $A$ in a different way. We will discuss this issue in Sec. III.4 where we solve the eigenvalue problem (4) on a finite grid by discretization. Here, we only remark that the indices $n$ and $l$ of Eq. (6) have the same meaning as in the probably more familiar case of a quantum mechanical particle in a spherical symmetric potential. Similar to the hydrogen atom, $n$ is related to the number of nodes of the collective mode, where $g(z)$ plays a role analogous to the wave function. Rotation symmetry implies a conserved angular momentum $l$ which enters the operator $A$ as a parameter. For $l=0$ the collective motion is isotropic, while for $l=1$ it is dipolar, and so on. Details can be found in Sec. III where we derive Eq. (5). We expect formula (6) to be applicable only for sufficiently small values of $n$ and $l$, because many nodes or a complicated angular structure may enter in conflict with the assumptions of hydrodynamics.

Apparently, for each $l>0$ there is a mode with $n=0$ and frequency $\sqrt{l}\omega_{0}$ which is independent of the equation of state. In particular, the lowest dipole mode is exactly at the trapping frequency. Up to now many measurements focused on the breathing and lowest quadrupole mode, $n=1,l=0$ and $n=1,l=2$, respectively. We emphasize that the measurement of two, three or more frequencies could achieve high precision when combined with theory: For a given equation of state an arbitrary number of eigenfrequencies can be obtained by our numerical procedure. Comparison to measurements in well-understood parameter regions tests the ability of a given theoretical method to compute the equation of state accurately. On the one hand, it will at least give us an estimate on the theoretical errors. On the other hand, if the agreement is very good, we can take advantage of a non-trivial dependence of the frequencies on thermodynamic quantities like temperature, density, chemical potential or equation of state parameters like $a$. For example, if we predict the lowest lying three monopole modes to have a specific temperature dependence, we may conclude from measuring these three modes in what temperature region we are.

Motivated by this we are looking for non-trivial relations between collective modes and thermodynamic quantities. We have already seen that a system purely described by mean field theory does not show such a behavior. We expect the mean field picture to be valid for very small gas parameter $na^{3}\ll 1$. However, as this parameter increases, higher order interaction effects become relevant in the equation of state. The leading order correction to the ground state energy density (2) at zero temperature has been calculated by Lee, Huang and Yang and is found to be LHY1

The corresponding pressure as a function of chemical potential to the same order in perturbation theory is given by

Apparently, the interaction strength can no longer be eliminated by a rescaling of $P(\mu)$ such that the ratio $P^{\mu}/P^{\mu\mu}$ entering the operator $A$ will always depend on $a$. Thus, the frequencies will depend on the coupling constant. While for a homogeneous system the thermodynamic observables are constant in space, for a confined gas one usually refers to their values at the center of the trap. Denoting the density in the center of the external potential by $n(0)$, $a^{3}n(0)$ provides a small parameter and one expects a shift of the breathing mode ($n=1,l=0$, i.e. $\omega_{B}=\sqrt{5}\omega_{0}$) PiSt2 ; Bra1

as compared to the Stringari formula above. We conclude that we may use the shift for small values of $n(0)a^{3}$ in order to determine either $n(0)$ or $a$ very precisely if the corresponding second quantity is known from another method.

From this example we can deduce a method to make high precision measurements using collective modes. If we assume the shifts of the frequencies $\omega_{n,l}$ to be continuous in $n(0)a^{3}$ or $\mu(0)a^{2}$, where $\mu(0)$ is the chemical potential in the center of the trap, then for small gas parameters, $n(0)a^{3}\ll 1$ or $\mu(0)a^{2}\ll 1$, the shift will be proportional to some power of the gas parameter. In the case of Eq. (9) we found this power to be $1/2$. However, this will depend on the system under consideration. Nevertheless, driving the system through a certain range of the gas parameter, e.g. by the use of a Feshbach Resonance, and then collecting the result in a double-logarithmic plot one will find the results to lie on a straight line for small gas parameter. This simple scaling behavior can be used for very precise measurements of $n(0)a^{3}$ after a proper calibration. The accuracy of this method is directly related to the accuracy in measuring the frequency shifts. In Fig. 1 we show the shifts beyond mean field due to the LHY-correction which we obtained by our numerical implementation described in Sec. III.4. We reproduce the prediction to a high precision in the regime of small gas parameter, where it is valid. Although this plot only shows the breathing mode we emphasize that it is possible to calculate the shifts of the whole frequency spectrum which are found to show a similar scaling behavior with $\sqrt{n(0)a^{3}}$. Fig. 1 shows that an experimental accuracy for frequency shifts better than $10^{-2}$ would be useful for exploring the regime of small gas parameters.

The calculation of higher order corrections to the LHY equation of state (7) is a difficult task. The next-to-next-to-leading order has been derived by Wu Wu1 ; Hua1 for hard-sphere bosons. The energy density receives two additional terms

Neglecting the logarithmic term, this leads to a correction which is proportional to $\mu^{3}a$ in Eq. (8). In Fig. 1 we have set $B=8$, see e.g. Ref. Nav1 for an overview over the expected values of $B$. Obviously, the general behavior of the frequency shifts is not influenced.

For large values of the gas parameter perturbation theory is no longer applicable and more sophisticated methods are necessary. The Functional Renormalization Group (FRG) for the effective average action Wet1 is a non-perturbative quantum field theoretical method and is in particular able to describe the full thermodynamics of systems which are realized in cold atom experiments. It does not rely on the expansion in a small parameter and therefore the most striking effects are expected in strongly interacting systems. For more information on the FRG see Refs. FRG . It is still under debate whether the regime $a\rightarrow\infty$ can be reached for a Bose gas at very low temperatures. One expects the condensate to get unstable then because of increasing importance of three-body losses.

In Fig. 1 we show the frequency shift obtained from the equation of state calculated with FRG at zero temperature Flo1 . We see that higher values of the gas parameter $n(0)a^{3}$ can no longer be identified with a unique $\delta\omega_{B}/\omega_{B}$. This is a fully non-perturbative effect. We recall that one of the assumptions in the beginning was that the interaction can be approximated to be point-like. However, this is only valid if the scattering length $a$ is much larger than the microscopic distance $\Lambda^{-1}$ where one can resolve the details of the interaction potential. For cold bosonic atoms this length is given by the typical range of the Van der Waals potential. An effective ultraviolet cutoff $\Lambda$ is necessary in any field theoretic treatment of the system where the contact interaction appearing in the Lagrangian of the non-relativistic system is not renormalizable in three dimensions. Therefore, the introduction of this scale is no relict of approximations but has a clear physical meaning. In perturbation theory one assumes $\Lambda$ to be infinity and thus the equation of state of the homogeneous system only involves one dimensionless parameter, the gas parameter $na^{3}$. A proper treatment has to account for the appearance of an additional length scale $\Lambda^{-1}$ such that two possible combinations, $na^{3}$ and $a\Lambda$, describe the system. As a result the shift of the breathing mode should rather be plotted as a two-dimensional surface depending on these two parameters. For small gas parameter the additional parameter $a\Lambda$ does not play a role – if we take equations of state $P(\mu,a)$ for different values of $a$ and then vary $n(0)$, all obtained shifts will lie on a line in the double-logarithmic plot. This is the expected scaling behavior. For higher densities it makes a difference whether we calculate the frequencies from $P(\mu,a,\Lambda_{1})$ or $P(\mu,a,\Lambda_{2})$ with $\Lambda_{1}\neq\Lambda_{2}$, as reflected by the spread of the points in Fig 1.

With regard to formula (6) we may wonder whether the cases $d=2$ and $d=1$ refer to highly anisotropic traps. These are of great relevance for experiments where one is often not dealing with a spherical symmetric potential. For example, the early experiments in the JILA group were performed in a disk-shaped confinement Jin1 while the MIT group used a cigar-shaped trap SK1 .

A really lower-dimensional system is obtained in a trap where the quantum gas is in its ground state in one or two directions. This will be the case for very tight confinement. When calculating the equation of state one can then neglect quantum and thermal fluctuations in these directions. The mean field Bose gas in this scenario is described by Eq. (6) when inserting $d=2$ or $d=1$ and $\alpha=1$. In these cases the system is isotropic in a lower-dimensional geometry.

A different situation arises if the system is not in its ground state in the directions of tight confinement and thus still three-dimensional or in a crossover between three and lower dimensions. Formula (6) cannot be applied in this case because the assumption of spherical symmetry is not justified. The solution of the hydrodynamic equations for these anisotropic cases is more involved.

We now consider the two-dimensional Bose gas. Its mean field equation of state is $\mu=g_{2D}n$ and the frequency of the breathing mode is found from Eq. (6) to be $\omega_{B}=2\omega_{0}$. The equation of state beyond mean field can be obtained by a Functional Renormalization Group approach Flo1 . The coupling constant $g_{2D}$, which satisfies $\mu=g_{2D}n$ for small $g_{2D}$, is dimensionless in two dimensions and shows a logarithmic running with the physical scale on which the experiment is performed. It will vanish for an infinitely large system. However, realistic experiments are always performed in traps so that in a harmonic potential the oscillator length $\ell_{0}=\sqrt{\hbar/m\omega_{xy}}$ provides the largest possible length scale of the physics under considerations. Therefore, when calculating the equation of state one only has to include quantum fluctuations with momenta bigger than $\ell_{0}^{-1}$, which acts as an infrared cutoff. This is also the reason why Bose–Einstein condensation is observed experimentally in two dimensions for small temperatures $0<T<T_{c}$. For an infinitely extended system the long-range order would be destroyed by fluctuations for all non-vanishing temperatures as required by the Mermin-Wagner theorem Mer1 .

In Figs. 2 and 3 we plot the shift of the breathing mode corresponding to $\omega_{B}=2\omega_{0}$ for the two-dimensional Bose gas at zero temperature with equation of state $P(\mu,g_{2D})$ from Functional Renormalization Group calculations Flo1 . Since the coupling constant is dimensionless there is no gas parameter for such a system. As we mentioned already in the three-dimensional case one always has a physical ultraviolet cutoff $\Lambda$ when dealing with a contact interacting in $d\geq 2$ dimensions. Since $\Lambda$ has dimension of inverse length, the dimensionless variable involving the chemical potential (or similar for the density) is $\mu(0)\Lambda^{-2}$, where $\mu(0)$ is the chemical potential of the gas in the center of the trap. A good choice for $\Lambda^{-1}$ is the oscillator length of tight trapping. We observe from Fig. 2 that the frequencies for small interactions depend only weakly on $\mu(0)\Lambda^{-2}$. For larger coupling $g_{2D}$ we find deviations from this behavior in Fig. 3. We emphasize that the dependence on the microscopic cutoff for two-dimensional bosonic systems is not a shortcoming of approximations, but a physical effect. It is due to the running of couplings. While the relevant coupling is dimensionless, the equation of state retains memory of scales in form of a dependence on the physical effective cutoff. This effect becomes substantial for large $g_{2D}$. Therefore, the strong dependence of the frequency shift on the dimensionless combination $\mu(0)\Lambda^{-2}$ for large interaction strength $g_{2D}\gtrsim 1$ in Fig. 3 signals the breakdown of universality. More microscopic properties beyond the single interaction parameter $g_{2D}$ are needed to describe a concrete experimental realization in this regime.

We arrive at two interesting experimental scenarios to be investigated. On the one hand, by measuring the collective modes one can distinguish whether one is working with a system which is still three-dimensional (disk, cigar) or a truly lower dimensional system. On the other hand, if the latter regime is reached experimentally it is of course tempting to verify the predictions of Fig. 3. The rather large negative frequency shifts can both test theoretical predictions like the occurrence of a minimum, and determine $\mu(0)$ or the corresponding density $n(0)$.

We now extend our considerations to non-vanishing temperature and allow for fluctuations of both density and temperature. If we stay below the critical temperature of Bose–Einstein condensation the system possesses both non-vanishing superfluid and normal fluid density and its hydrodynamic behavior has to be described by two-fluid hydrodynamics Lan1 ; Kha1 ; Put1 .

As explained in Sec. III, we again find an eigenvalue problem which has to be solved in order to get the collective frequencies of the system. It has the general form

where the differential operators $A$, $B$, $C$ and $D$ are defined in Eqs. (73) - (79). These operators depend on the equation of state $P(\mu,T)$ and the normal fluid density $n_{n}(\mu,T)$. Both quantities have to be provided by an underlying microscopic theory. The structure of the eigenvalue problem is similar to the simple zero temperature case. We will show below that the zero temperature limit of Eq. (11) is indeed given by Eq. (4) and this behavior will also be found in the frequencies.

We see that calculating the frequencies at non-vanishing temperature is as straightforward as it was at $T=0$. However, results on the temperature dependence of hydrodynamic collective modes are rare in the literature Ho1 . This is because the equation of state and normal fluid density are in general not known or only in certain ranges of temperature and chemical potential. But when applying the local density approximation $\mu\rightarrow\mu-V_{ext}(r)$ one drives through a wide interval of values for the chemical potential and therefore a complete equation of state in $\mu$ has to be provided. Our motivation is in applying the Functional Renormalization Group, which is capable of providing the full thermodynamics of cold quantum gases, especially the functions $P(\mu,T)$ and $n_{n}(\mu,T)$. For the three-dimensional weakly interacting Bose gas the equation of state has been calculated by Lee and Yang LY1 . Their formulas are valid for small gas parameter $na^{3}$ for temperatures not too close to the critical temperature. We comment on the critical behavior later in this section. The Lee–Yang equation of state is summarized in App. B. We neglect the next-to-leading order LHY-correction (7) because we are interested in thermal effects.

In Ho1 Shenoy and Ho have calculated the temperature behavior of collective modes of a trapped Bose gas from two-fluid hydrodynamics and Lee–Yang theory for $0.6<T/T_{c}<1.2$. The authors considered the range of parameters where in addition to $na^{3}\ll 1$ also $an\lambda_{T}^{2}\ll 1$ holds. ($\lambda_{T}=(2\pi\hbar^{2}/mk_{B}T)^{1/2}$ is the thermal wavelength.) In this limit the integral expressions for the equation of state from Lee–Yang theory can be evaluated analytically. We will, however, keep the full expressions for the equation of state.

Before discussing the collective modes of a Bose gas at non-vanishing temperature we comment on some aspects of the static configuration of the trapped gas which are necessary for the interpretation of our results. For a homogeneous system with temperature $T$ we can calculate the critical chemical potential $\mu_{c}(T)$ of the superfluid phase transition. Within the local density approximation in the trapped system there will be a radius $r_{c}$ corresponding to the phase boundary between the superfluid and the normal regions of the cloud. This radius fulfills $\mu_{c}=\mu(0)-V_{ext}(r_{c})$ with the chemical potential $\mu(0)$ in the center of the trap. The characteristic picture of the density profile consists in a narrow condensate peak in the center of the trap which is surrounded by a broad thermal cloud of the normal gas. Of course, there is also a non-vanishing contribution of the normal component to the inner regions. We visualize this behavior in Fig 4.

It is now apparent that for the description of a trapped gas at any nonzero temperature one has to know the equation of state for both the superfluid and normal phase. In particular the presence and dynamics of the normal gas are of importance. There can be oscillations of the thermal cloud itself and, furthermore, it provides a non-trivial background for the oscillations of the condensate. As we approach zero temperature the thermal cloud vanishes. Condensate oscillations correspond to solutions $\delta n$ of the hydrodynamic equations which are only nonzero inside a sphere with radius $r_{c}$.

The critical temperature $T_{c}$ of the trapped gas is defined as the temperature where the condensate peak appears in the center of the cloud. This can be reformulated as $\mu(0)=\mu_{c}$ or equivalently $r_{c}=0$. For Lee–Yang theory the critical temperature is given by

If we set $\mu(0)=\mu_{c}$ and calculate the number of particles $N$ in the trap, we find the ideal gas prediction $T_{c}^{(0)}(N)=\hbar\omega_{0}(N/\zeta(3))^{1/3}$ to be satisfied for small gas parameter. There are several possibilities to change the parameter $T/T_{c}$. For fixed $\mu(0)$, the particle number will change when decreasing the temperature. Another way of going to low temperatures is to keep the particle number fixed. Then the chemical potential in the center of the trap has to be adjusted appropriately. If non-destructive measurements of frequencies become possible, a particular promising experimental way would be a change of $T/T_{c}$ by a smooth variation of the trap frequency, keeping the number of trapped atoms fixed.

In Fig. 5 the isotropic ($l=0$) temperature dependent oscillation frequencies in an isotropic harmonic trap are shown for fixed particle number $N=10^{6}$ and $a/\ell_{0}=0.0005$, where $\ell_{0}=\sqrt{\hbar/m\omega_{0}}$ is the oscillator length. This particular choice of parameters implies a small gas parameter $na^{3}$, which is required in order to use Lee–Yang theory, but does not guarantee the applicability of hydrodynamics. Nevertheless, we expect this approach to capture the characteristics of the temperature dependence of hydrodynamic collective modes. Due to thermodynamic constraints the equation of state for the strongly coupled case will share common features with Lee–Yang theory and thus will show a qualitatively similar oscillation spectrum. For a quantitatively precise description of the collective modes of the weakly coupled system, which is not the purpose of this paper, one needs to include corrections to the hydrodynamic behavior by means of collisionless dynamics. Our results then constitute a limiting case.

We observe a rich spectrum of frequencies. The oscillations can be classified as condensate oscillations and oscillations of the thermal cloud. The condensate oscillations correspond to the branches which disappear at $T_{c}$. The computability of frequencies above the critical temperature by our method is a manifestation of the fact that the Landau two-fluid model remains formally valid for vanishing superfluid component. Of course, the frequencies for $T\geq T_{c}$ could be computed equivalently with a one-fluid model for the thermal liquid. We find the oscillation frequencies of the thermal cloud to be very close to the spectrum which one obtains using the equation of state for an ideal Bose gas and vanishing superfluid density. In particular, the lowest mode which also exists above $T_{c}$ is given by $\omega=2\omega_{0}$, see for example Ref. Str2 .

For low temperatures ($T/T_{c}\lesssim 0.1$) the two-fluid region, defined by $r\leq r_{c}$ with $r_{c}$ from above, starts to grow considerably and the normal fluid in the outer region vanishes. However, we were not able to obtain the eigenfrequencies of this crossover numerically in a sufficiently accurate manner. In order to get an understanding of the very low temperature limit of the collective modes we therefore used an implementation where the cloud is cut off at $r_{c}$ and thus only the two-fluid region is left. One then obtains a spectrum of the type shown in Fig. 8, which shows a smooth limit for $T\rightarrow 0$, as is discussed below. (For increasing $T$ there remains a systematic error in Fig. 8 because with increasing temperature the normal mantle surrounding the two-fluid core of the cloud cannot be neglected anymore.)

From Fig. 5 it is clear that there are level-crossing frequencies over wide ranges of $T/T_{c}$. These branches can be used as a thermometer by measuring a certain number of frequencies and then locating them in the plot. How this can be done experimentally by response techniques is explained in Sec. II.4.

An apparent feature of Fig. 5 is the appearance of bifurcations points where one branch of collective excitation frequencies splits up into two distinct branches. The meaning of these points is very intuitive. While for low temperatures superfluid and normal fluid component are closely connected, the modes split up at intermediate values of $T/T_{c}$. One of them represents the superfluid oscillation, which vanishes at the critical temperature, the other one corresponds to the oscillation of the normal fluid which remains present also above $T_{c}$. The black solid dots in Fig. 5 correspond to points before and after a splitting. We indeed find that for low temperatures there is only one solution $\delta\mu$ and $\delta T$ for this particular eigenvalue, which then smoothly splits up into a set of two solutions with different properties. We make this statement clear in Figs. 6 and 7 where the (isotropic) oscillations of chemical potential and temperature are shown for the frequencies corresponding to the black dots in Fig 5, distinguished by the labels “upper” and “lower” branch, meaning higher and lower frequency, respectively.

One has to keep in mind, that the oscillations of particle number density and entropy density, $\delta n$ and $\delta s$, vanish at infinity because the thermodynamic functions in Eqs. (56) and (57) vanish. However, the oscillations of density and entropy density are related to $\delta\mu$ and $\delta T$ via

where the matrix contains the higher derivatives of the equation of state $P(\mu,T)$ applied to local density approximation. (We use here the notation of Sec. III.) Since this matrix vanishes at infinity, $\delta\mu$ and $\delta T$ can be nonzero. At zero temperature, where we have a sharp cloud radius of the condensate, the fluctuations of the chemical potential $\delta\mu$ are discontinuously going to zero at this radius.

From our discussion of the equilibrium density profile we see that for all temperatures below the critical temperature $T_{c}$ there is a radius $r_{c}$ where the critical equation of state $P(\mu_{c},T)$ has to be known. The vicinity of the point $\mu_{c}$ might be small but the behavior of the thermodynamic functions in this interval could nevertheless influence the frequencies significantly. Indeed, the superfluid phase transition is of second order and we expect the specific heat at constant pressure $C_{P}=T(\partial S/\partial T)_{P,N}$ and the isothermal compressibility $\kappa_{T}=-V^{-1}(\partial V/\partial P)_{T,N}$ to be singular at $\mu_{c}$. Their behavior in the critical region as a function of temperature is $C_{P}\sim|T-T_{c}|^{-\alpha}$ and $\kappa_{T}\sim|T-T_{c}|^{-\gamma}$ with the critical indices PV of the three-dimensional XY-universality class. The specific heat shows a cusp while the isothermal compressibility diverges. A systematic study of these effects on the collective modes would be very interesting. Lee–Yang theory does not allow for such an investigation.

Using the Lee-Yang formulas given in App. B we can calculate numerically the coefficient functions appearing in the operators $A$, $B$, $C$, and $D$ in Eq. (11). For very low temperatures we observe that $B$ and $C$ vanish while $A$ approaches its zero temperature expression, given by Eq. (5),

This corresponds to a decoupling of superfluid and thermal degrees of freedom. This is expected to happen in the zero temperature limit since the few atoms in the thermal cloud have no influence on the dynamics of the condensate. The operator entering the eigenvalue problem is now block diagonal and all zero temperature frequencies are reproduced. They correspond to Stringari’s mean field formula given above. (Note that we have neglect the LHY-correction.) It is interesting that $D$ does not vanish. Indeed, the coefficient functions $a_{1}$, $d_{1}$ and $d_{2}$ entering $A$ and $D$ from Eqs. (77) and (79) satisfy

This can be seen by evaluating the integrals for $P(\mu,T)$ and $n_{n}(\mu,T)$ numerically or applying an expansion for $T\ll\mu$ using only the phonon part of the spectrum, see App. D. It can be shown that $a_{1}/d_{1}$ corresponds to the ratio of the squared first and second velocities of sound, $c_{1}^{2}$ and $c_{2}^{2}$, respectively. Therefore we find the well-known relation

to be satisfied for very low temperatures.

We show in App. D that the behavior of the coefficient functions of $D$ as $T\rightarrow 0$, Eq. (15), can be used to derive the eigenvalues of $D$ by slight modifications from the eigenvalues of $A$, which are given by Eq. (6) with $\alpha=1$. The corresponding eigenfrequencies are found to be

The experimental verification of this formula will be difficult, however. Although for low $T$ (and $T\rightarrow 0$) these frequencies are solutions to the two-fluid hydrodynamic equations, there are only very few atoms left in the normal phase to oscillate. From Eq. (17) it is apparent that for $n=1,l=0$ one obtains $\sqrt{2/3}\omega_{0}$ which is below the trapping frequency. This cannot be achieved by formula (6) with an effective polytropic index $\alpha_{eff}$. The existence of such modes below the trapping frequency might be a special feature of superfluid hydrodynamics just as the appearance of a second velocity of sound. Indeed, in their early experiments Stamper-Kurn et al. SK1 used a cigar-shaped harmonic potential with trapping frequency $\omega_{0z}/2\pi=18.04(1)$ Hz in axial direction. They observed an out-of-phase oscillation between condensate and thermal cloud of a Bose gas at $T=1\mu K$ in the hydrodynamic regime. The critical temperature was around $1.7\mu K$. The measured frequency of axial motion was found to be $\omega_{0}/2\pi=17.26(9)$ Hz, which is below $\omega_{0z}$.

In Fig. 8 we show our numerical results for the frequencies in the zero temperature limit using the full Lee–Yang equation of state for $\mu(0)/\hbar\omega_{0}=10$ fixed and $a/\ell_{0}=0.0005$. We find formulas (6) and (17) to be in perfect agreement with our data.

We have seen that information is contained both in the precise values of single modes and in the whole spectrum of the oscillating system. While the lowest lying modes can be excited by varying the parameters of the trapping potential (e.g. the potential depth) and then measuring the frequency in a free expansion, the whole spectrum could be obtained by an other method. In principle one can fit a superposition of arbitrary many modes to the freely expanding cloud but the contribution of the higher modes will only be weak. We suggest to use a response measurement instead to locate the frequencies of the trapped cloud. In order to keep the notation short we restrict this discussion to the zero temperature case.

The eigenfrequencies and corresponding eigenfunctions obtained so far represent normal coordinates of an oscillating system, the trapped cold gas. The equation of motion $\partial_{t}^{2}\delta\mu+A\delta\mu=0$ has solutions $\delta\mu_{nlm}=\bar{g}_{nl}(r)r^{l}Y_{lm}e^{-i\omega_{nl}t}$ which are harmonic oscillations with frequency $\omega_{nl}$. (More precisely, one should use the operator $E$ defined below Eq. (122) instead of $A$. This is of no relevance here, because they have the same eigenvalues.)

To describe an experimental situation we need to include dissipation effects into our description. We therefore introduce a phenomenological damping term by replacing $\partial^{2}_{t}\rightarrow\partial^{2}_{t}+\Gamma\partial_{t}$, with damping constant $\Gamma$. The new eigenfrequencies of the system will be denoted by $\Omega_{nl}$. We use Fourier decomposition and write $\delta\mu(t)=\int\mbox{d}\Omega\delta\mu(\Omega)e^{-i\Omega t}$. Since we already know the eigenvalues of $A$ we immediately obtain the quadratic equation

for possible solutions $\Omega_{nl}$. The latter simplify due to the fact that $t_{2}=1/\Gamma$ represents a characteristic time scale over which damping takes place. In order to observe oscillations this scale should be substantially larger than $t_{1}=1/\omega_{nl}$. The other choice would correspond to the overdamped case. Hence $\omega_{nl}^{2}\gg\Gamma^{2}$ and we have

for the new eigenfrequencies of the system. Due to damping effects they acquire an imaginary part and their real part is slightly shifted. This phenomenological treatment gives no explanation of the microscopic nature of the damping terms. It would be interesting to derive the relation between the damping constants $\Gamma$ and the five transport coefficients $\eta,\kappa_{T},\zeta_{1},\zeta_{2}$ and $\zeta_{3}$ of two-fluid hydrodynamics Put1 . In principle, one expects $\Gamma$ to depend on the mode in the form of some continuous function $\Gamma(\Omega)$. For simplicity, we use here the same phenomenological $\Gamma$ for all modes.

The picture of the trapped gas as an oscillating system will provide us with an intuition how individual (and thus also higher lying) modes can be addressed by virtue of response techniques. Consider a classical point particle with Hamiltonian $H=\frac{1}{2m}p^{2}+\frac{m}{2}\omega^{2}x^{2}$. With a damping term included its equation of motion will have the form of Eq. (18) and is therefore independent of the mass of the particle. However, when applying a small external force $f(t)$ the equation of motion for the Fourier components of $x(t)$ becomes

Obviously, the response $\chi(\Omega)=x(\Omega)/f(\Omega)$ is proportional to the inverse of the mass $m$. It is therefore more difficult to excite a heavy particle than a lighter one. We introduce an analog quantity $m_{nlm}$ which represents the “mass” for the collective oscillation $\delta\mu_{nlm}$. We call $m_{nlm}$ the “response coefficients”. The relation between the amplitude of the oscillations of the atom cloud and the oscillation of the external perturbation $f(\Omega)$ is governed by $m_{nlm}$.

Let us consider a periodic perturbation of the harmonic trapping potential,

The perturbation is assumed to be small in comparison to $V_{ext}$. As we show in App. E this small variation acts as a driving force on the oscillating system (18). There it is also shown that in order to excite an oscillation with angular dependence $\delta\mu_{nlm}\propto Y_{lm}$, $\delta V_{ext}$ has to have the same angular dependence. In Eq. (132) we give a general formula for the response coefficient $m_{nlm}$ of the collective mode. Given $\Omega\simeq\omega_{nl}$ we then have

for the imaginary part of the response function defined in Eq. (131). The latter is related to dissipation of energy in the system.

For a demonstration we restrict ourselves to the case of a three-dimensional mean-field Bose gas, i.e. $P(\mu)\propto\mu^{2}$, and consider only the isotropic modes with $l=0$. In addition we assume a monomial perturbation $F(r)=r^{q}$ with $q\geq 0$. This includes the cases of a spatially constant background field $\delta V_{ext}(\vec{x},t)=\mbox{cos}(\Omega t)$ ($q=0$) and a time-variation of the trapping frequency, $\delta V_{ext}=r^{2}\mbox{cos}(\Omega t)$ ($q=2$). Interestingly, we find in App. E that not all values of $q$ can be used to excite the whole spectrum. For example, there is no response at all to $q=0$, while for $q=2$ only the breathing mode will be excited and for $q=4$ only the lowest two modes oscillate. We show this behavior in Fig. 9. (In our intuitive picture developed above the absence of response corresponds to an infinite mass.) The normalization in Fig. 9 is such that the lowest mode always has the same amplitude. The relative amplitudes, both for oscillations with different frequencies, and for oscillations with a given frequency, but for different excitations, can be used as a further observable for testing the system. The computation of ratios of amplitudes for different $\Omega$ needs knowledge of $\Gamma(\Omega)$. In principle, this can be measured from the width of the resonance. Ratios of amplitudes for fixed $\Omega$ and different radial or angular dependence of the perturbation are independent of $\Gamma(\Omega)$ and $m_{nlm}$. In principle, amplitudes as in Fig. 9 can be measured both by dissipation ($\mbox{Im}\chi$) or by the amplitude of the cloud oscillations ($\mbox{Re}\chi$).

We conclude that the variation of the external potential with a certain frequency and radial and angular dependence allows for an excitation of individual modes. The oscillation can then be detected by imaging the density profile or, hopefully, by new non-invasive methods.

We mentioned in the beginning of this section that we expect the features of the Bose gas to be generic for other systems of cold atoms. At nonzero temperature this requires the applicability of the two-fluid model which is related to a two-component order parameter. The latter can be written as $\phi(\vec{x})=|\phi(\vec{x})|e^{i\theta(\vec{x})}$ at each space-time point. If $\theta(\vec{x})$ is varying slowly in space one can define the superfluid velocity as its gradient, $\vec{v}_{s}\propto\nabla\theta$. Thus our description remains valid for systems which are described by a complex scalar field which acquires a non-vanishing expectation value. This includes in particular composite bosons like the atom-atom dimers on the BEC-side of the BEC-BCS crossover. For an imbalanced system it would then be interesting to include the conserved spin degrees of freedom of the remaining fermions as an additional macroscopic variable in the hydrodynamic equations.

At zero temperature any system with hydrodynamic equations of motion might be considered. An example can be found in Ref. Cit1 where the collective oscillation frequencies of a one-dimensional dipolar quantum gas have been calculated.

The polytropic equation of state $P\propto\mu^{\alpha+1}$ is of particular interest for scale-invariant systems like ultracold Fermi gases at the unitary point where the scattering length $a$ diverges. At zero temperature only the chemical potential remains as a dimensionful quantity, the pressure has to be proportional to $\mu^{(d+2)/2}$. Indeed, introducing the Bertsch parameter $\xi$ we can write the chemical potential at the unitary point as $\mu=\xi\varepsilon_{F}$, where $\varepsilon_{F}\propto n^{2/d}$ is the Fermi energy density. Using $\mbox{d}P=n\mbox{d}\mu$ we arrive at $P(\mu)\propto\mu^{(d+2)/2}$ as claimed. This time we know the equation of state exactly, but only for zero temperature. At the unitary point and for $T=0$ the oscillation frequencies obey the exact Eq. (6), with $\alpha=3/2$ for a three-dimensional spherical trap.

Around the unitary point and in the dilute, three-dimensional regime, the equation of state can be parametrized as Bul1

with an additional parameter $\zeta$ implying a shift

of the breathing mode.

Observations on quantum fluids suggest that there is a critical temperature $T_{c}$ such that for $0\leq T\leq T_{c}$ the macroscopic hydrodynamic motion can be divided into two kind of flows, a superfluid (subscript $s$) and a normal fluid (subscript $n$) one. Their velocity fields $\vec{v}_{s}(\vec{x},t)$ and $\vec{v}_{n}(\vec{x},t)$, respectively, enter the momentum density $\vec{g}=n_{s}\vec{v}_{s}+n_{n}\vec{v}_{n}$ with coefficients satisfying $n=n_{s}+n_{n}$, where $n(\vec{x},t)$ is the density of the gas or fluid under consideration. We might be lead to the conclusion that the substance is built out of two different kinds of fluids. However, the true nature of superfluidity consists in quantum physics and our description is only formally correct. The two flows can be distinguished by the following properties. The superfluid velocity field is a gradient field and hence irrotational. Furthermore, entropy is only carried by the normal fluid part. The full hydrodynamic equations under these constraints were set up by Lev. D. Landau in a famous and still very worth reading paper in 1941 Lan1 . We concentrate on ideal hydrodynamics where dissipation terms are neglected and entropy is conserved. In the homogeneous case, the underlying conservation laws read

with stress tensor $\Pi_{ik}=n_{n}v_{n,i}v_{n,k}+n_{s}v_{s,i}v_{s,k}+P\delta_{ik}$, entropy density $s$, pressure $P=-\varepsilon+Ts+\mu n+n_{n}(\vec{v}_{n}-\vec{v}_{s})^{2}$, chemical potential $\mu$ and internal energy density in the rest frame of the gas, $\varepsilon$. These equations have to be supplemented by the equation of state $P(\mu,T)$ and the normal fluid density $n_{n}(\mu,T)$. For dissipative hydrodynamics one would replace Eq. (27) by the conservation of energy and include an additional term in Eq. (28). We use dimensionless units described in appendix A, setting $\hbar=k_{B}=m=1$.

In order to understand what changes when applying an external field, we first derive the local density approximation from a thermodynamic point of view. For this purpose consider a gas contained in a potential $V_{ext}(\vec{x})$ such that local equilibrium is reached at each point of space. If we take two neighboring small but yet macroscopic parts of the gas their energies $E_{1}$ and $E_{2}$ and particle numbers $N_{1}$ and $N_{2}$ will adjust in a manner maximizing the entropy $S=S_{1}+S_{2}$ under the constraints $E_{1}+E_{2}=E$ and $N_{1}+N_{2}=N$, respectively. This implies temperature and the full chemical potential to be constant inside the trap. However, the full chemical potential is given by the Gibbs free energy per particle and thus reads $\mu_{full}(\vec{x})=\mu_{hom}(P(\vec{x}),T)+V_{ext}(\vec{x})$, where $\mu_{hom}(P,T)$ is the equilibrium chemical potential of the homogeneous system as a function of pressure $P$ and temperature $T$. We conclude that a system where the chemical potential $\mu$ in the homogeneous equilibrium functions is substituted according to $P_{hom}(\mu,T)\rightarrow P_{hom}(\mu-V_{ext},T)$ etc. behaves like a system trapped in an external potential of large spatial extend. This rule is called local density approximation.

An external potential can be included in Eqs. (26) and (28) by virtue of force terms,

The static solution denoted by a subscript $0$ gives the expected expressions $\mu_{0}(\vec{x})+V_{ext}(\vec{x})=\bar{\mu}$ and $\nabla P_{0}(\vec{x})+n_{0}(\vec{x})\nabla V_{ext}(\vec{x})=\vec{0}$ of static equilibrium. Using the Gibbs-Duhem relation $\mbox{d}P=s\mbox{d}T+n\mbox{d}\mu$ and thus $\nabla P_{0}=s_{0}\nabla T_{0}+n_{0}\nabla\mu_{0}$ at every space point, we conclude $s_{0}\nabla T_{0}=\vec{0}$ in equilibrium. We recognize the right thermodynamic behavior to emerge naturally. Pressure and densities of mass, entropy and energy are constant in time but space dependent, while temperature is constant all over the trap and the chemical potential follows the rule of local density approximation. The parameter $\bar{\mu}$ can be adjusted in order to get a certain particle number $N$. The fluid velocities $\vec{v}_{s}$ and $\vec{v}_{n}$ vanish in equilibrium.

The next step towards our eigenvalue problem consists in expanding the two-fluid equations in small fluctuations around their static equilibrium solution such that time-dependent quantities only appear linearly. We write

In the same way we linearize $P=P_{0}+\delta P,n=n_{0}+\delta n,s=s_{0}+\delta s$. Eqs. (25), (27), (29) and (30) then become

with $\delta\vec{g}=n_{s,0}\delta\vec{v}_{s}+n_{n,0}\delta\vec{v}_{n}$. Using the relation $\nabla P_{0}=s_{0}\nabla T_{0}+n_{0}\nabla\mu_{0}$ we get $\nabla\delta P=s_{0}\nabla\delta T+\delta n\nabla\mu_{0}+n_{0}\nabla\delta\mu=s_{0}\nabla\delta T-\delta n\nabla V_{ext}+n_{0}\nabla\delta\mu$, since $\nabla T_{0}=\vec{0}$ and $\nabla\mu_{0}=-\nabla V_{ext}$. This can be used to cast Eq. (36) into the form

In the following we will consider a spherically symmetric harmonic trapping potential $V_{ext}(\vec{r})=\frac{m}{2}\omega_{0}^{2}r^{2}$. The generalization of our formalism to other radial potentials $V_{ext}(\vec{r})=V_{ext}(r)$ is straightforward. In Eq. (69) we give a formulation of the eigenvalue problem which holds for arbitrary external potentials $V_{\rm ext}(\vec{x})$ and thus also allows for a treatment of anisotropic trapping geometries.

We assume the fluctuations of the physical quantities to be a harmonic oscillation in time with frequency $\omega$,

and similar expressions for the other time-dependent functions. The remaining spatial deviations are again denoted by $\delta\mu$ etc. Since time has completely dropped out of our system of partial differential equations, this slight abuse of notation hopefully does not lead to confusion. We recognize our obtained set of equations

to be an eigenvalue problem for a differential operator acting on $\delta\mu$, $\delta T$, $\delta\vec{v}_{s}$ and $\delta\vec{v}_{n}-\delta\vec{v}_{s}$.

We derive the zero temperature eigenvalue problem in terms of an equation of state given in the form $P(\mu)$. We introduce the notation $\mu_{0}(r)=\mu(0)-V_{ext}(r)=\mu(0)-\omega_{0}^{2}r^{2}/2$. Since we have set $\hbar=k_{B}=m=1$, we can express all dimensionful quantities with respect to $\omega_{0}$, see App. A. In the following we set $\omega_{0}=1$ in all formulas, which means that we are measuring time in units of the inverse trapping frequency. Thus, $\mu_{0}(r)=\mu(0)-r^{2}/2$. Furthermore, we define

where $\mu(0)$ is the chemical potential in the center of the trap. We can write $\delta n(\vec{x})=P^{\mu\mu}_{0}(r)\delta\mu(\vec{x})$ because $\mbox{d}P=n\mbox{d}\mu$. There are no temperature fluctuations present at $T=0$ and thus the only degrees of freedom are $\delta\mu$ and $\delta\vec{v}_{s}$ described by Eqs. (40) and (41). The latter equations can be decoupled. Writing $n_{0}(r)=P^{\mu}_{0}(r)$ we arrive at the Stringari wave equation Str1

So far, $\delta\mu$ depends on the $d$-dimensional spatial coordinate $\vec{x}$. Since we are dealing with a spherically symmetric trap we can classify the possible solutions to the eigenvalue problem via an ansatz

where $f_{lm}$ are spherical harmonics. (For a definition of spherical harmonics in $d\leq 3$ dimensions we refer to App. C. Note that we have $l=0,1$ for $d=1$.) When applying this ansatz to Eq. (45) we benefit from the facts that $\nabla\mu_{0}\cdot\nabla f_{lm}=0$ and $\partial_{r}P(\mu_{0}(r))=-rP_{0}^{\mu}(r)$. In addition we use the relations