# The Surface of a Bose-Einstein Condensed Atomic Cloud

## I Introduction

Properties of trapped clouds of Bose-Einstein condensed atoms have been investigated intensively both experimentally and theoretically over the past few years, since the experimental realization of Bose-Einstein condensation in dilute atomic gases. At zero temperature, the behavior of the order parameter, $`\psi `$, is determined by the time-dependent Gross-Pitaevskii equation, a Schrödinger equation with the nonlinear term proportional to $`|\psi |^2\psi `$ added:

$$\frac{\mathrm{}^2}{2m}^2\psi (𝐫,t)+V(𝐫)\psi (𝐫,t)+U_0|\psi (𝐫,t)|^2\psi (𝐫,t)=i\mathrm{}\frac{\psi (𝐫,t)}{t}.$$
(1)   
Here $`U_0`$ is the effective two-particle interaction, which may be expressed in terms of the scattering length $`a`$ according to $`U_0=4\pi a\mathrm{}^2/m`$, where $`m`$ is the mass of an atom. The confining potential is denoted by $`V(𝐫)`$. The static structure of the cloud is determined by the time-independent Gross-Pitaevskii equation, which reads

$$\frac{\mathrm{}^2}{2m}^2\psi (𝐫)+V(𝐫)\psi (𝐫)+U_0|\psi (𝐫)|^2\psi (𝐫)=\mu \psi (𝐫),$$
(2)   
where $`\mu `$ is the chemical potential. For clouds containing a sufficiently large number of atoms with repulsive interactions, many properties may be calculated to a good approximation using the Thomas-Fermi approximation, in which one neglects the kinetic energy term in the Gross-Pitaevskii equation. One then finds a density profile

$$|\psi (𝐫)|^2=(\mu V(𝐫))/U_0,\mu V(𝐫),|\psi (𝐫)|^2=0,\mu <V(𝐫).$$
(3)   
The Thomas-Fermi approach may also be used for calculating collective modes, and this has been done by Stringari.

A number of phenomena are associated with the surface region of the cloud. One is that the contribution of the kinetic energy term to the total energy of the cloud comes mainly from the surface region. Also Stringari has identified surface modes of oscillation in the Thomas-Fermi approximation. To understand such phenomena in detail, it is useful to consider the properties of a planar surface, and to approximate the trapping potential by a linear function of the coordinates. The potential is then given by

$$V(𝐫)=Fx,$$
(4)   
where the coordinate $`x`$ measures distances in the direction of $`V`$. The presence of the $`^2`$ term in the Gross-Pitaevskii equation leads to a rounding-off of the density profile in the surface region over a distance of order

$$\delta =\left(\frac{\mathrm{}^2}{2mF}\right)^{1/3},$$
(5)   
and within this approach the static structure of the surface region has been studied, and contributions to the total energy calculated .

The purpose of this paper is to consider properties of plane surfaces of Bose-Einstein condensed clouds in a linear potential. This is of interest both for giving analytical expressions for a number of properties of large clouds, as well as for giving physical insight. First, we consider the static properties of the surface region, and show how the kinetic energy term in the Gross-Pitaevskii equation gives rise to an effective surface tension. We also consider the density profile in the vicinity of the surface, and show that the Gross-Pitaevskii equation and the Thomas-Fermi approximation to it lead to identical results for the column density of atoms down to a point well within the cloud. In the second part of the paper we consider surface collective modes. We begin by considering the Thomas-Fermi approach, which yields a mode with the dispersion relation $`\omega ^2=(F/m)q`$, where $`\omega `$ is the angular frequency and $`q`$ the magnitude of the wave vector for the mode. This is exactly the same form as for gravity waves on the surface of a fluid, with the role of gravity being played by the trapping potential. At higher wave numbers there are contributions to the mode frequency arising from the $`^2\psi `$ term in the Gross-Pitaevskii equation, and we evaluate the effect of these by employing a trial wave function that allows for the rounding-off of the order parameter profile in the surface region. The extra contributions to $`\omega ^2`$ are proportional to $`q^4\mathrm{ln}(1/q\delta )`$, and they will be shown to have a ready interpretation in terms of the effective surface tension introduced in the discussion of static properties.

## II Static properties

In this section we introduce the concept of an effective surface tension, and consider the number of particles that are associated with the rounding-off of the density profile in the surface region. We shall take the origin of the $`x`$-coordinate to be at the point where the potential is equal to the chemical potential, which is where the density vanishes in the Thomas-Fermi approximation.

Let us now calculate the contribution to the energy coming from the kinetic energy term in the Hamiltonian. As was discussed in Ref. , the total kinetic energy may be written in terms of either $`_{\mathrm{}}^{\mathrm{}}(\psi ^{})^2𝑑x`$ or $`_{\mathrm{}}^{\mathrm{}}\psi \psi ^{\prime \prime }𝑑x`$, where the prime denotes differentiation with respect to $`x`$: the kinetic energy density is not a well defined quantity. For calculating the total energy either expression may be used, but it must be used consistently. In the case of the surface problem under investigation, we wish to be able to associate a contribution to the kinetic energy with the surface, and there is no unambiguous way of doing this, since one could use either of the energy expressions integrated to some point well inside the surface (rather than $`\mathrm{}`$ as in the case of the total kinetic energy). However, these differ by a constant. For definiteness, we shall employ the symmetrized expression and define the kinetic energy per unit area of the surface as

$$ϵ_K=\frac{\mathrm{}^2}{2m}_L^{\mathrm{}}(\psi ^{})^2𝑑x,$$
(6)   
where the point $`x=L`$ is chosen to lie well inside the surface. From Ref.  this may be seen to be given by

$$ϵ_K\frac{\mathrm{}^2}{8m}\frac{F}{U_0}\mathrm{ln}\left(\frac{L}{0.240\delta }\right),$$
(7)   
where the coefficient 0.240 was found by numerical integration. If we had chosen the other expression for the kinetic energy density, the result would differ by $`\mathrm{}^2F/4mU_0`$. The coefficient of $`\mathrm{ln}L`$ is independent of the choice of kinetic energy density, but terms independent of $`L`$ are not. This difference, however, is unimportant in application to surface modes, where the logarithmic term is the one of most interest. If this expression were independent of $`L`$, the kinetic energy would have precisely the same form as a surface tension, since the energy would be proportional to the area of the surface. Because the kinetic energy density falls off only slowly away from the surface, the total kinetic energy increases logarithmically with $`L`$, but this dependence is sufficiently mild that it is still useful to think of the kinetic energy as being physically analogous to a surface tension.

Next we consider the density profile, and compare the Thomas-Fermi result with the exact one. For $`|x|\delta `$ the Thomas-Fermi wave function is a good approximation to the exact one and thus $`\psi ^{\prime \prime }`$ is negative. Consequently, as one can see from the Gross-Pitaevskii equation, the density is depressed below the Thomas-Fermi value. Outside the cloud, the density is increased, due to the quantum-mechanical tail of the wave function. We shall now demonstrate the somewhat surprising result that the number of particles associated with the depression of densities inside the cloud is exactly equal to the number associated with regions of increased density.

To prove this result it is simplest to consider the momentum density, $`𝐠`$, defined by

$$g_i=\frac{\mathrm{}}{2i}(\psi ^{}\frac{\psi }{x_i}\psi \frac{\psi ^{}}{x_i}),$$
(8)   
whose time dependence may be determined from the Gross-Pitaevskii equation (1). By differentiating (8) and using (1) we find that the momentum density satisfies the conservation condition

$$\frac{g_i}{t}+\frac{T_{ij}}{x_j}=n\frac{V}{x_i},$$
(9)   
where $`n(=|\psi |^2)`$ is the density. Here the stress tensor $`T_{ij}`$ is given by

$$T_{ij}=T_{ij}^0+P\delta _{ij},$$
(10)   
where

$$T_{ij}^0=\frac{\mathrm{}^2}{2m}\left(\frac{\psi }{x_i}\frac{\psi ^{}}{x_j}+\frac{\psi ^{}}{x_i}\frac{\psi }{x_j}\frac{1}{2}\frac{^2|\psi |^2}{x_ix_j}\right)$$
(11)   
is the free-particle stress tensor, while $`P`$ is the pressure arising from the interaction,

$$P=\frac{1}{2}n^2U_0.$$
(12)   
For the one-dimensional problem we are considering, the condition for hydrostatic equilibrium is simply

$$\frac{T_{xx}}{x}=Fn,$$
(13)   
from which it follows that

$$T_{xx}(L)=F_L^{\mathrm{}}n𝑑x,$$
(14)   
since the stress tensor vanishes for large positive $`x`$. Deep within the cloud the order parameter varies as $`(x)^{1/2}`$, and therefore $`T_{xx}^0(L)1/L`$, which vanishes for $`L\mathrm{}`$. Consequently

$$_L^{\mathrm{}}n𝑑x=\frac{1}{2}n^2(L)U_0+𝒪(1/L).$$
(15)   
The first term is the result one obtains in the Thomas-Fermi approximation, where $`n=xF/U_0`$ for $`x<0`$, and thus one sees that the total column density outside a point well within the surface is the same as that in the Thomas-Fermi approximation, apart from corrections of order $`1/L`$. Physically the result follows from the condition that the stress tensor at any point must be balanced by the force on all the material at larger values of $`x`$, a result familiar in the context of equilibrium of fluids and gases in the presence of gravitational fields.

## III Surface modes of oscillation at long wavelengths

We turn now to time-dependent situations and consider oscillations. For this purpose it is convenient to use instead of (1) an equivalent set of equations for the density, and the local velocity of the condensate. The equations of motion are obtained by writing $`\psi `$ in terms of its amplitude $`f`$ and phase $`\varphi `$, $`\psi =fe^{i\varphi }`$. The number density is given by $`n=f^2`$, while the velocity $`𝐯`$ is given by $`𝐯=\mathrm{}\varphi /m`$. By inserting $`\psi =fe^{i\varphi }`$ into (1) and separating the equation into real and imaginary parts one obtains the two equations

$$\frac{(f^2)}{t}=\frac{\mathrm{}}{m}(f^2\varphi ),$$
(16)   
which is the equation of continuity,

$$\frac{n}{t}+(n𝐯)=0,$$
(17)   
and

$$\mathrm{}\frac{\varphi }{t}=\frac{\mathrm{}^2}{2mf}^2f+\frac{1}{2}mv^2+V(𝐫)+U_0f^2.$$
(18)   
We eliminate the phase variable by taking the gradient of (18), using $`𝐯=\mathrm{}\varphi /m`$. The resulting equation is written as

$$m\frac{𝐯}{t}=(\delta \mu +\frac{1}{2}mv^2),$$
(19)   
where

$$\delta \mu =V+U_0n\frac{\mathrm{}^2}{2m\sqrt{n}}^2\sqrt{n}\mu _0.$$
(20)   
Since it is the gradient of $`\delta \mu `$, which enters the acceleration equation (19), we are free to subtract a constant from $`\delta \mu `$. We have in (20) chosen to subtract the value of the equilibrium chemical potential $`\mu _0`$, which implies that $`\delta \mu `$ is zero in equilibrium, i. e. under stationary conditions. The equation $`\delta \mu =0`$ is the time-independent Gross-Pitaevskii equation.

In the Thomas-Fermi approximation, the kinetic energy term is neglected, and the equilibrium density is therefore given by

$$n_0U_0+V(x,y,z)=\mu _0.$$
(21)   
Within the Thomas-Fermi approximation one also neglects the kinetic energy term involving $`\delta n`$ in the expression (20) for $`\delta \mu `$. This yields

$$\delta \mu =U_0\delta n.$$
(22)   
The Thomas-Fermi approximation for the modes should be a good one provided the wavelength of the mode is large compared with the healing distance $`\delta `$. With these approximations we may readily linearize the equations (17) and (19), and eliminate $`\delta \mu `$ by means of (22). The result is

$$m\frac{^2\delta n}{t^2}=U_0(n_0\delta n).$$
(23)   
If we only consider oscillations with the time dependence $`\delta n\mathrm{exp}(i\omega t)`$, the differential equation (23) simplifies to

$$\omega ^2\delta n=\frac{U_0}{m}(n_0\delta n+n_0^2\delta n).$$
(24)   
We investigate the surface modes in a two-dimensional configuration with the linear ramp potential considered in the previous section. In the $`y`$\- and $`z`$-directions there is translational invariance, and therefore the solution must have the form of plane waves for these coordinates. We denote the wavenumber of the mode by $`q`$, and take the direction of propagation to be the $`z`$-axis. In the Thomas-Fermi approximation the condensate density in equilibrium, $`n_0`$, is then given by $`n_0(x)=Fx/U_0`$ for $`x<0`$, while it vanishes for $`x>0`$. It follows that (24) has a solution of the form

$$\delta n=Ae^{qxiqz},$$
(25)   
which describes a wave propagating on the surface, and decaying exponentially in the interior. Since (25) satisfies $`^2\delta n=0`$, while the gradient of the equilibrium density is given by $`(F/U_0,0,0)`$, we obtain by inserting (25) into (24) the dispersion relation

$$\omega ^2=\frac{F}{m}q.$$
(26)   
This has the same form as for a gravity wave propagating on the surface of an incompressible ideal fluid in the presence of a gravitational field $`g=F/m`$.

The solution (25) is however not the only one which decays exponentially in the interior. To investigate the solutions to (24) more generally we insert a function of the form

$$\delta n=f(qx)e^{qx+iqz},$$
(27)   
and obtain the following second order differential equation for $`f(y)`$,

$$y\frac{d^2f}{dy^2}+(2y+1)\frac{df}{dy}+(1ϵ)f=0,$$
(28)   
where $`ϵ=\omega ^2/gq`$. By introducing the new variable $`z=2y`$ one sees that Eq. (28) becomes the differential equation for the Laguerre polynomials $`L_n(z)`$, provided $`ϵ1=2n`$. We have thus obtained the general dispersion relation for the surface modes

$$\omega ^2=\frac{F}{m}q(1+2n),n=0,1,2\mathrm{}$$
(29)   
with the associated density oscillations given by

$$\delta n(x,z,t)=AL_n(2qx)e^{qx+iqzi\omega t},$$
(30)   
with $`A`$ being an arbitrary constant.

To make contact with Stringari’s calculation, we note that he found the dispersion relation of modes in an isotropic harmonic trap to be given by $`\omega ^2=\omega _0^2[l(1+2n)+3n+2n^2]`$, where $`l`$ is the angular momentum quantum number and $`n`$ the radial one, which gives the number of nodes in the radial direction. For $`l`$ much greater than 1, the dispersion relation becomes $`\omega ^2=\omega _0^2l(1+2n)`$. The wavenumber of the mode at the surface of the cloud is given by $`q=l/R`$, and therefore the dispersion relation is $`\omega ^2=\omega _0^2qR(1+2n)`$, which is precisely the same as the result (29) we obtained above, since the force due to the trap at the surface of the cloud is simply $`F=\omega _0^2R`$ per unit mass. For large values of $`l`$ it is thus a good approximation to replace the harmonic oscillator potential by the linear ramp, as one might expect since the surface modes are concentrated within a distance of order $`R/l`$ from the surface. It should be noted that the $`n=0`$ mode frequencies for the plane surface with a linear ramp potential agree with the frequencies of the nodeless radial modes (corresponding to $`n=0`$) for a harmonic trap at all values of $`l`$. For modes with radial nodes ($`n0`$), the two results agree only when $`l`$ is much greater than $`n`$.

## IV Surface modes at shorter wavelengths

When the healing length $`\delta `$ is not negligible compared with the wavelength, there are corrections to the dispersion relation. In the case of gravity waves on the surface of a liquid, modes at shorter wavelengths are affected by the surface tension, and the dispersion relation is given by $`\omega ^2=gq+\sigma q^3/\rho `$, where $`\sigma `$ is the surface tension, and $`\rho `$ is the mass density of the fluid. We now explore modes on the surface of a Bose-Einstein condensed cloud at shorter wavelengths, and we shall show that there are contributions to $`\omega ^2`$ of order $`q^4\mathrm{ln}(1/q\delta )`$ which may be understood in terms of the effective surface tension introduced in the discussion of static properties.

The basic problem is to solve Eq. (19) including the quantum pressure term in the expression for the chemical potential. Rather than attacking the problem directly, which leads to two coupled second-order differential equations, we shall adopt a variational approach, which will allow us to calculate the leading corrections to the Thomas-Fermi result for the mode frequencies for small $`q`$.

In order to determine the dispersion relation of surface modes at shorter wavelengths we employ a trial wave function that allows us to calculate the total energy in terms of two variables which describe the displacement of the surface and the local velocity, respectively. In terms of these variables the energy functional assumes the form of that of a harmonic oscillator, from which we may extract the frequency as a function of $`q`$.

The trial wave function is motivated by the solution found above in the Thomas-Fermi approximation. To lowest order we may describe the motion of the surface in a traveling wave by modifying the ground state wave function in two respects. First one shifts the spatial variable, thereby allowing for displacements of the surface, and one introduces a phase factor to take into account motion of the particles. Explicitly, the wave function is given by

$$\psi (x,z,t)=\psi _{TF}(x\mathrm{\Delta }(x,z,t))\mathrm{exp}i\varphi \psi _{TF}(x)\mathrm{\Delta }(x,z,t)\psi _{TF}^{}(x)+i\varphi \psi _{TF}(x).$$
(31)   
Here $`\mathrm{\Delta }(x,z,t)=\xi _0\mathrm{exp}qx\mathrm{cos}(qz\omega t)`$, with $`\xi _0`$ constant, and $`\varphi =(mv_0/\mathrm{}q)e^{qx}\mathrm{sin}(qz\omega t)`$, $`v_0`$ being the amplitude of the velocity at the surface of the cloud, given by

$$𝐯=\frac{\mathrm{}\varphi }{m}=v_0(\mathrm{sin}(qz\omega t),0,\mathrm{cos}(qz\omega t)).$$
(32)   
Let us now turn to the more general case. We expect that at frequencies small compared with the characteristic frequency $`\mathrm{}/2m\delta ^2`$ associated with adjustments of the density profile in the region within $`\delta `$ of the surface, the density profile will be able to adjust essentially instantaneously to its equilibrium form corresponding to the local number of particles per unit area, even if the Thomas-Fermi approximation is not valid. The real part of the wave function in the vicinity of the surface is thus of the equilibrium form, but with a possible translation perpendicular to the surface. We shall therefore use a trial function which has the same form as for the Thomas-Fermi case, but with the equilibrium Thomas-Fermi wave function replaced by the exact one. From this we shall calculate the energy, and evaluate oscillation frequencies. We write the wave function in the form

$$\psi =\psi _0+\delta \psi ,$$
(33)   
where $`\delta \psi `$ is the part due to the oscillation. For the present purposes it is simplest to consider a standing wave, and therefore we adopt the following form

$$\delta \psi (x,t)=(\xi (t)\psi _{0}^{}{}_{}{}^{}(x)+i\psi _0(x)\varphi _0(t))e^{qx}cos(qz).$$
(34)   
The $`x`$-component of the velocity of the surface is given either in terms of the time derivative of the surface displacement, or in terms of the $`x`$-derivative of the phase of the wave function. This leads to the consistency condition $`\dot{\xi }=\mathrm{}q\varphi _0/m`$.

We use the trial function (34) to evaluate the energy functional

$$E=𝑑𝐫\left[\frac{\mathrm{}^2}{2m}|\psi |^2+V(𝐫)|\psi (𝐫,t)|^2+\frac{1}{2}U_0|\psi (𝐫,t)|^4\right],$$
(35)   
to second order in $`\xi `$ and $`\dot{\xi }`$. The zeroth order term gives rise to an unimportant constant, while the first order term vanishes, since the trial function is a solution of the Gross-Pitaevskii equation. The interesting physics is contained in the second order term $`E^{(2)}`$. The $`\xi ^2`$-part of the kinetic-energy contribution contains an integral over $`x`$ of the form

$$_{\mathrm{}}^{\mathrm{}}𝑑x(\psi _{0}^{}{}_{}{}^{}e^{qx})^{}(\psi _{0}^{}{}_{}{}^{}e^{qx})^{}=_{\mathrm{}}^{\mathrm{}}𝑑xe^{2qx}\psi _0^{}(\psi _{0}^{}{}_{}{}^{\prime \prime \prime }+2q\psi _{0}^{}{}_{}{}^{\prime \prime }+q^2\psi _{0}^{}{}_{}{}^{}).$$
(36)   
The third-order derivative in (36) is eliminated by differentiating the Gross-Pitaevskii equation,

$$\frac{\mathrm{}^2}{2m}\psi _0^{\prime \prime \prime }=V^{}\psi _0V\psi _{0}^{}{}_{}{}^{}3U_0\psi _0^2\psi _{0}^{}{}_{}{}^{}.$$
(37)   
The last two terms on the right hand side of (37) yield contributions that are canceled by those coming from the potential and interaction energies in (35). The remaining terms may be combined using $`V^{}=F`$ and partial integration, and yield for the energy per unit area

$$E^{(2)}=\frac{\xi ^2}{4}\left[Fq_{\mathrm{}}^{\mathrm{}}𝑑xe^{2qx}\psi _{0}^{}{}_{}{}^{2}+\frac{\mathrm{}^2}{m}_{\mathrm{}}^{\mathrm{}}𝑑x\left(q^2e^{2qx}\psi _0^2\right)\right]+\frac{m\dot{\xi }^2}{4}_{\mathrm{}}^{\mathrm{}}𝑑xe^{2qx}\psi _{0}^{}{}_{}{}^{2}.$$
(38)   
This is of the same form as for an harmonic oscillator, $`E^{(2)}=\frac{1}{2}C_1\xi ^2+\frac{1}{2}C_2\dot{\xi }^2`$ where $`C_1`$ and $`C_2`$ are constants, and the frequency is given by $`\omega ^2=C_1/C_2`$. Integrating by parts the term involving $`F`$, we obtain the final result

$$\omega ^2=\frac{F}{m}q+\frac{\mathrm{}^2q^4}{m^2}I(q),$$
(39)   
where the dimensionless quantity $`I(q)`$ is given by

$$I(q)=\frac{_{\mathrm{}}^{\mathrm{}}𝑑xe^{2qx}(\psi _{0}^{}{}_{}{}^{})^2}{q^2_{\mathrm{}}^{\mathrm{}}𝑑xe^{2qx}\psi _{0}^{}{}_{}{}^{2}}.$$
(40)   
The first term in (39) gives the frequency of the surface mode in the Thomas-Fermi approximation, while the second term involving $`I(q)`$ is a correction term.

We may evaluate the leading long-wavelength corrections to the dispersion relation $`I(q)`$ for $`q1/\delta `$ by splitting up the range of integration into two regions

$$_{\mathrm{}}^{\mathrm{}}𝑑xe^{2qx}(\psi _{0}^{}{}_{}{}^{})^2=_{\mathrm{}}^L𝑑xe^{2qx}(\psi _{0}^{}{}_{}{}^{})^2+_L^{\mathrm{}}𝑑xe^{2qx}(\psi _{0}^{}{}_{}{}^{})^2,$$
(41)   
where $`\delta L1/q`$. In the first of the two integrals we may use the Thomas-Fermi approximation $`\psi _0=\sqrt{Fx/U_0}`$, since $`L\delta `$. In the second we may replace $`\mathrm{exp}2qx`$ by 1, since $`ql1`$. By using partial integration one finds

$$_L^{\mathrm{}}𝑑x(\psi _{0}^{}{}_{}{}^{})^2=\frac{F}{2U_0}_L^{\mathrm{}}𝑑x\psi _{0}^{}{}_{}{}^{\prime \prime }\psi _0$$
(42)   
since the value of $`\psi _0\psi _0^{}`$ at $`x=L`$ may be evaluated in the Thomas-Fermi approximation, resulting in $`F/2U_0`$. The integral on the right hand side of (42) occurs in calculations of the kinetic energy associated with the surface region, and it has been evaluated previously. Its asymptotic form is $`(F/4U_0)\mathrm{ln}(L/1.776\delta )`$, which corresponds to the kinetic energy contribution in Thomas-Fermi theory cut off at a distance $`1.776\delta `$, where the coefficient was determined by numerical integration. We obtain consequently

$$_{\mathrm{}}^l𝑑xe^{2qx}(\psi _{0}^{}{}_{}{}^{})^2=\frac{F}{4U_0}(2\mathrm{ln}2qL\mathrm{ln}\gamma )+\frac{F}{4U_0}\mathrm{ln}\frac{L}{1.776\delta },$$
(43)   
where $`\gamma (1.778)`$ is the Euler constant. The denominator in (40) is evaluated in the Thomas-Fermi approximation, resulting in $`F/4U_0`$. We finally obtain

$$I(q)\mathrm{ln}q\delta +\mathrm{ln}(e^2/3.552\gamma )\mathrm{ln}q\delta +0.15.$$
(44)   
Thus we find

$$\omega ^2\frac{F}{m}q+\frac{\mathrm{}^2q^4}{m^2}\left[\mathrm{ln}q\delta +0.15\right].$$
(45)   
The qualitative behavior is easy to understand by analogy with the surface tension contribution to the frequency of gravity waves on the surface of a fluid, where the contribution to $`\omega ^2`$ is $`\sigma q^3/\rho `$. In the present problem, the effective surface tension depends logarithmically on the length scale, which is given by the wavelength. The $`q^4`$-dependence exhibited by (45) then results from dividing $`q^3`$ by the effective density in the region where the fluid is moving, this being of order the fluid density at the distance $`1/q`$ from the surface, or $`F/U_0q`$. Thus one sees that the surface tension is very weakly dependent on the trap parameters and the atomic scattering length, which occur only in the logarithm.

This result agrees to the order indicated with the result of using the full solution $`\psi _0(x)`$ to the (equilibrium) Gross-Pitaevskii equation in evaluating the expression (40). In a recent paper, Fetter and Feder analyzed the corrections to the excitation frequencies in the Thomas-Fermi limit for an atomic cloud confined by a spherically symmetric trap. By employing the matching conditions of boundary-layer theory, they were able to demonstrate that the leading correction to the Thomas-Fermi limit is of order $`R^4`$, where $`R`$ is the Thomas-Fermi radius, while terms of order $`R^4\mathrm{ln}R`$ were found to be absent. Repeating their analysis for the different geometry which we are considering, we have explicitly verified that the matching conditions allow for the presence of terms of order $`q^4\mathrm{ln}q`$, which we have found in the present paper.

This result for the surface mode frequency can be obtained in a more rigorous fashion by a variational approach. The resulting equations of motion of $`\xi `$ and $`\varphi _0`$ take the form of that of a classical harmonic oscillator with frequency given by (39). The details of this calculation are described in Appendix A.

We thank Emil Lundh for helpful contributions.

## A variational approach

The Gross-Pitaevskii equation may be derived from the variational principle

$$\delta 𝑑tL=0,$$
(A1)   
where

$$L=𝑑𝐫\frac{i\mathrm{}}{2}\left(\psi ^{}\frac{\psi }{t}\psi \frac{\psi ^{}}{t}\right)E.$$
(A2)   
Here $`E`$ is the energy functional given by equation (35). For $`\psi `$ we adopt our ansatz, (33), and the equations of motion for $`\xi `$ and $`\varphi _0`$ may then be determined from the variational principle. To second order in $`\xi `$ and $`\varphi _0`$ the Lagrangian is

$$L=\frac{\mathrm{}}{4}I_1(\xi \dot{\varphi }_0\varphi _0\dot{\xi })(E^{(0)}+E^{(2)}),$$
(A3)   
where $`E^{(0)}`$ is the ground-state energy obtained by inserting $`\psi =\psi _0`$ in equation (35), while $`E^{(2)}`$ is given by

$$E^{(2)}=\frac{\xi ^2}{4}\left[FqI_2+\frac{\mathrm{}^2}{m}I_3\right]+\frac{\mathrm{}^2q^2}{4m}\varphi _0^2I_2,$$
(A4)   
in terms of the integrals

$$I_1=\frac{1}{2}_{\mathrm{}}^{\mathrm{}}𝑑xe^{2qx}(\psi _0^2)^{^{}},$$
(A5)   
$$I_2=_{\mathrm{}}^{\mathrm{}}𝑑xe^{2qx}\psi _0^2,$$
(A6)   
and

$$I_3=_{\mathrm{}}^{\mathrm{}}𝑑x(q^2e^{2qx}\psi _{0}^{^{}}{}_{}{}^{2}).$$
(A7)   
The resulting equations of motion for $`\xi `$ and $`\varphi _0`$ are

$$\frac{\mathrm{}}{2}I_1\dot{\varphi }_0\frac{1}{2}\xi (FqI_2+\frac{\mathrm{}^2}{m}I_3)=0$$
(A8)   
and

$$\frac{\mathrm{}}{2}I_1\dot{\xi }+\frac{\mathrm{}^2q^2}{2m}I_2\varphi _0=0.$$
(A9)   
By partial integration of (A5) we see that $`I_1=qI_2`$, and subsequently, by eliminating $`\varphi _0`$ between the last two equations, one obtains the result (39).