{'fraction_non_alphanumeric': 0.0380291865337379, 'fraction_numerical': 0.013755237682415836, 'mean_word_length': 4.119230769230769, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Vortex Dynamics Near the Surface of a Bose-Einstein Condensate

The center-of-mass dynamics of a vortex in the surface region of a Bose-Einstein condensate is investigated both analytically using a variational calculation and numerically by solving the time-dependent Gross-Pitaevskii equation. We find, in agreement with previous works, that away from the Thomas-Fermi surface, the vortex moves parallel to the surface of the condensate with a constant velocity. We obtain an expression for this velocity in terms of the distance of the vortex core from the Thomas-Fermi surface that fits accurately with the numerical results. We find also that, coupled to its motion parallel to the surface, the vortex oscillates along the direction normal to the surface around a minimum point of an effective potential.


I. INTRODUCTION
Following the experimental observation of vortices in Bose-Einstein condensates [1,2,3,4,5,6,7,8,9], vortex dynamics has been a subject of intensive interest [10,11,12,13,14,15,16,17,18,19,20,21,22]. Theoretical studies have been performed on the center-of-mass motion of the vortex core in the bulk of the condensate [10,15,22]. References [10,22] have studied the dynamics of the core taking into account the dissipation from the thermal cloud. Dissipation causes the vortex to spiral out of the condensate giving it a finite lifetime in agreement with the experimental observations [2]. Once they reach the surface of the condensate, the vortices decay into surface excitations. It is also believed that vortices enter the condensate from the surface at the point when surface excitations become unstable [19,21,23,24,25,26,27]. Hence, it is worth investigating the details of the dynamics of a vortex in the surface region of the condensate.
E. Lund et al. [28] have pointed out that near the surface of the condensate, the quadratic trapping potential can be approximated by a potential that is linear in distances normal to the surface. As a result of this approximation, the equilibrium density profile will be also linear in distances normal to the surface of the condensate except for distances very close to the Thomas-Fermi surface. An immediate question then arises: How will the vortex move in such an inhomogeneous density background? The fact that the density is linear makes it possible to perform analytical treatment to the problem provided that the size of the core of the vortex is much less than the distance over which the background density varies significantly. This question has been addressed extensively by Anglin [20], where the author used a boundary-layer theory to calculate both the density profile and phase of the vortex, from which he obtained the equilibrium velocity of the vortex. It turned out that the component of the velocity in the direction normal to the surface of the condensate vanishes, while the component parallel to the surface is a constant that depends on the distance between the core of the vortex and the surface. In this paper we address basically the same problem but using a variational calculation. The main aim is to demonstrate that the variational calculation performed here generates the previous knowledge about the vortex motion near the surface and captures the main features of it. We go one step ahead of previous works by obtaining a rather more accurate expression for the component of the vortex velocity parallel to the surface. Furthermore, we show that, within few coherence lengths from the Thomas-Fermi surface, the vortex oscillates along the direction normal to the surface around a minimum point of an effective potential. These oscillations are coupled to its motion parallel to the surface and disappear when the vortex is away from the surface.
In our variational calculation we use an ansatz wavefunction that takes into account the equilibrium density profile and phase of both the background and the vortex. The variational parameters are the coordinates of the vortex core and the center-of-mass velocity components normal to the surface of the condensate and parallel to it. The equations of motion for these parameters are then derived from a lagrangian that corresponds to the time-dependent Gross-Pitaevskii equation.
Using a linear density profile for the background is an approximation that breaks down at the Thomas-Fermi surface. Our variational calculation is, thus, expected to be accurate only away from the Thomas-Fermi surface. Furthermore, we neglect the internal dynamics of the core of the vortex. This is also justified away from the Thomas-Fermi surface where the size of the core of the vortex is almost constant.
In the following, we start, in subsection II A, by writing down the ansatz wavefunction. Using this wavefunction, we calculate the lagrangian from which we derive the equilibrium properties in subsection II B and the equations of motion for the variational parameters in subsection II C. In Section III, we present the results of the numerical simulation of the vortex motion, and compare them with those of the variational calculation. In Section IV we summarize our main conclusions.

A. The Lagrangian
The dynamics of the time-dependent Bose-Einstein condensate is described by the Gross-Pitaevskii equation Here g is the effective two-particle interaction which is proportional to the s-wave scattering length a according to g = 4πah 2 /m, where m is the mass of an atom, and µ is the equilibrium chemical potential. The harmonic trapping potential V (r ′ ) is given by where ω 0 is the characteristic frequency of the trap. Here, we take for simplicity V (r ′ ) to be isotropic. For distances near the surface of the condensate, the quadratic trapping potential V (r ′ ) can be approximated by a linear potential [28,29] where x is a coordinate normal to the surface of the condensate, R is the Thomas-Fermi radius given by V (R) = µ, and F = mω 2 0 R is a force constant. The last equation is obtained by using the transformation r ′ = r + Rx in Eq. (2). This suggests a planner geometry in which the surface of the condensate is approximated by an infinite plane. The x-coordinate will be normal to the plane while the y and z coordinates will be parallel to it. The origin of this coordinate system is then shifted from the center of the condensate to the surface. As a result of this shift, the bulk of the condensate occupies the region x < 0 and the surface region will be near x = 0. This is illustrated in Fig. 1.
In this planner geometry, the Gross-Pitaevskii equation takes the form The equilibrium density profile n 0 (x) = |ψ(x)| 2 is obtained from the last equation by setting the time derivative to zero For large negative values of x, the kinetic energy operator can be neglected [28,30] leading to the Thomas-Fermi equilibrium density n TF 0 (x) = −(F/g)x. We consider a vortex with a core located at r = x 0 (t)x + y 0 (t)ŷ and axis along the z-direction as shown in Fig. 1. The vortex has a center-of-mass velocityv = v x (t)x + v y (t)ŷ. Our trial wavefunction is given by Here, the first factor on the right hand side is n TF 0 (x). The second factor corresponds to the square root of the density of a vortex in a uniform background. The function Φ(ρ, φ) in the exponential is the phase of the vortex that takes into account the inhomogeneouty of the background density. The last factor represents the center-of-mass velocity of the vortex. The polar coordinates ρ = (x − x 0 (t)) 2 + (y − y 0 (t)) 2 and φ = tan −1 ((x − x 0 )/(y − y 0 )) have the core of the vortex as their origin, as indicated in Fig. 1. Near the surface, the Thomas-Fermi approximation breaks down and the density background is nonlinear. The kinetic energy starts to become important leading to an exponential decay of the density. Our trial wavefunction does not take this nonlinearity into account. Therefore, the results of this variational calculation are not expected to be accurate in the surface region where the Thomas-Fermi approximation breaks down.
For the phase of the vortex we use an expression similar to the one derived by Anglin [20] The first term of this expression represents the usual φ-dependance of a vortex in a uniform background. The second term results from the linear inhomogeneouty of the background. It should be noted that the angle φ defined in this paper is the complement of that defined in Ref. [20]. This expression is obtained using a boundary-layer theory which is described here briefly. Vortices have structure in their density over a small length scale of the order of the coherence length ξ 0 = 1/ 8πan 0 (x 0 ), and phase that extends over the whole system. Here, ξ 0 is the coherence length calculated with the background density at x 0 , and δ = (h/mF ) 1/3 is the surface depth [28]. In the boundary-layer theory, the vortex problem is first solved exactly in the inner region ρ ∼ ξ 0 and the outer region ρ ≫ ξ 0 . The inner asymptotics of the outer solution and the outer asymptotics of the inner solution are then obtained by expanding the outer and the inner solutions in powers of two perturbation parameters characteristic to the two regions. Finally, the two solutions are matched to determine the solution in the intermediate region.
The velocity and phase of the vortex are obtained as a result of this matching procedure. The phase, given by Eq. (7), corresponds to the inner asymptote of the outer solution. Thus, this expression is valid only near the vortex and is not expected to be accurate for ρ ≫ ξ 0 . However, in the present calculation, we take this expression to be valid for all values of ρ. We account for the fact that it is not valid for distances away from the vortex core, by introducing in the argument of the logarithm the parameter α, which works as an effective cutoff on the logarithmic term. The value of α is then determined by fitting the velocity calculated variationally in subsection II B with the one obtained from the numerical solution of the time-dependent Gross-Pitaevskii equation in section III. The vortex density profile is taken as the square of [31] χ(ρ) = ρ It turns out in subsection II C that the equations of motion of y 0 (t) and v y (t) are independent of χ. Therefore, specifying the functional form of χ is not necessary for the y-component of the motion. Due to the nonlocal nature of the vortex excitation, the energy of the vortex diverges logarithmically with the size of the system [32]. To be able to perform the variational calculation for a finite and constant number of atoms N , we consider in calculating the lagrangian only the atoms within a cylinder of radius b and length l such that the cylinder axis coincides with the vortex axis. The fact that we restrict the number of atoms within the cylinder to be constant, is guaranteed by requiring ψ(r, t) to be normalized to N . The normalized wavefunction takes the form where The lagrangian that corresponds to Eq. (4) is given by where ∆E[ψ, ψ * ] is the energy associated with the presence of the vortex. This energy is obtained by subtracting the energy of the vortex-free background from the vortex energy [32], namely where and the vortex-free energy is The terms on the right hand side of the last equation represent the mean-field and trapping potential energies of the vortex-free background, respectively. These were calculated using the wavefunction The prefactor of ψ vf guarantees its normalization to N . For the subtraction of the background energy from the vortex energy to be correct, it is essential that the number of atoms of both systems to be the same.
Using the trial wavefunction Eq. (9) in the lagrangian of Eq. (11), we show in Appendix A that the lagrangian per atom takes the form where the coefficients N 1 , . . . , N 7 are functions of b, and N 2 and N 4 have also x 0 dependence ( See Appendix A.). The dot on y 0 , v x , and v y denotes derivative with respect to time. The parameter γ = gN/N 1 l is an energy that characterizes the average mean-field energy since N/N 1 l represents the average density within the cylinderical volume surrounding the vortex. This lagrangian will be the basis of the calculations of the equilibrium and nonequilibrium properties in the next two subsections.

B. Equilibrium Properties
The equilibrium values of the variational parameters can be obtained by minimizing the energy functional with respect to the variational parameters. The energy functional can be readily obtained from Eq. (16) by setting time derivatives to zero and then multiplying by -1, namely Minimizing ∆E with respect to v x gives v x = 0. As we shall see in the next subsection, v x is proportional toẋ 0 in the limit x 0 ≪ −δ. Therefore, x 0 will be a constant throughout the motion of the vortex. In other words, the vortex will be moving parallel to the surface of the condensate in agreement with previous works [16,20] and the results of the numerical simulation of section III. Minimizing the energy with respect to v y yields In the hydrodynamic limit, b ≫ ξ, the coefficient The quantityh/2mx 0 is the equilibrium velocity of a vortex in a homogeneous density background located at a distance x 0 from a hard wall and is moving parallel to it. The logarithmic term corresponds, therefore, to the contribution of the linear inhomogeneouty in the background density. The value of αb is determined by fitting this expression to the numerical values of v y (t). This is performed in section III where it turns out that αb ≈ 0.11. Thus Eq. (19) reads This is to be compared with the expression derived by Anglin [20] using a boundary-layer theory, namely v y (t) = (h/2mx 0 (t)) (1.96 + (3/2) ln (x 0 /δ)). In Fig. 2, we plot v y as a function of x 0 using Eq. (20) (20) represents one of the main results of this paper as it shows that the variational calculation does indeed lead to a rather accurate description of the vortex velocity parallel to the surface of the condensate.
Substituting the equilibrium expression of v y from Eq. (18) in Eq. (17), we obtain an energy functional that is a function of v x and x 0 , namely where is a function of x 0 that can be considered as an effective potential for the vortex. Here, we have neglected a term containing N 2 4 since it is one order of magnitude smaller than the logarithmic term of Eq. (7). Inspection shows that this effective potential has a minimum for negative values of x 0 . This means that the vortex can acquire oscillations around this minimum. Taking this into consideration and the fact that for large x 0 , we have v x =ẋ 0 , suggests following a collective coordinates approach to find the frequency of oscillation [29]. In this approach the energy functional, Eq. (21), is put in the form of that of a simple harmonic oscillator The frequency of oscillation will be simply ω = k/m. This is achieved by expanding the coefficients N 1 , N 2 (x 0 ), N 5 , N 6 , and N 7 in powers of ξ/b keeping terms up to the second order, namely: In this limit, the effective potential reduces to which has a minimum at Expanding U (x 0 ) around x eq , we obtain The energy functional becomes The frequency of vortex oscillation is then given by Noting that γ =h 2 /2mξ 2 ∼ n(x 0 ) and that in the limit |x 0 | ≫ δ, the Thomas-Fermi approximation gives n(x 0 ) ∼ |x 0 |, we conclude that ω ∼ |x 0 | −4 ln b/ξ. Thus, for large x 0 , these vortex oscillations disappear and x 0 will be constant as we have also found above by minimizing ∆E with respect to v x . In the next subsection, we derive equations of motion that describe the vortex motion and oscillations.

C. Equations of Motion
The equations of motion follow from the Euler-Lagrange equations of the lagrangian Eq.
The prime on N 2 (x 0 ) denotes a derivative with respect to x 0 . In obtaining the last equation in terms of x 0 only, we have used Eqs. (29) and (30). In addition, we have neglected terms containing N 2 4 (x 0 ) and N 4 (x 0 )N ′ 4 (x 0 ) since these terms are of an order of magnitude less than the logarithmic term in the phase Φ(ρ, φ) of Eq. (7).
From the first two equations, we obtain quite generallÿ This is a general result in the sense that it does not depend on the functional form of χ, the upper limit of the radial integration b, the cutoff of the logarithmic part of the phase α, or x 0 . The coordinates of the vortex are then described completely by Eqs. (31) and (32).
In the limit |x 0 | ≫ δ, Eqs. (30) and (31) reduce tov x = x 0 andẍ 0 = 0, respectively. Therefore, x 0 will be in general a linear function of t, namely x 0 (t) = c 1 t + c 2 , where c 1 and c 2 are constants. We have shown in the previous subsection that v x = 0 for the energy to be minimum. This results in that c 1 = 0, and thus x 0 is a constant as we have also found in the previous subsection. We notice also that, in this limit, Eq. (29) reduces to v y =ẏ 0 . This implies that in the limit |x 0 | ≫ δ, the variational parameters v x and v y are the conjugate velocities of x 0 and y 0 , respectively.
In the limit b ≫ ξ, the coefficients N 1 , . . . , N 7 can be expanded in powers of ξ/b, as was shown in the previous subsection, and Eq. (31) simplifies to where U (x 0 ) is the effective potential of the vortex defined in Eq. (23). Expanding the last equation around x eq and using Eq. (25), we get Using the initial conditions x 0 (0) = x eq + X 0 andẋ 0 (0) = 0, the solution of the last equation reads where Ω = ω/(1 − b 2 /4x 2 eq ) and ω is given by Eq. (27). Notice that in the limit b ≫ δ, the quantity b 2 /x 2 eq → 0, which leads to Ω → ω. Equation (31) can be solved numerically for given values of b. We plot this numerical solution in Fig. 3 together with analytical solution Eq. (35). This figure shows that the analytical solution is a reasonable approximation for the exact solution of Eq. (31).
To ensure that these oscillations are not an artifact of our variational calculation, we have solved numerically the time-dependent Gross-Pitaevskii equation for initial vortex distances from the Thomas-Fermis surface that are of order δ. We found that the vortex indeed oscillates in the x-direction as predicted by the variational calculation. Furthermore, we show next that we can obtain good agreement between the variational calculation and the numerical one for both the frequency and amplitude of the oscillations. In the present variational calculation, the amplitude of vortex oscillation X 0 and frequency ω are essentially determined by Eqs. (24) and (27). Scaling, length to δ and energy toh 2 /mδ 2 = F δ, these two equations reduce, in the zeroth order of ξ/b, to: x eq = (δ/ξ) 2 δ and ω/(F/mδ) = (ξ/δ) 4 ln (b/ξ), where we have usedh 2 /2mξ 2 = γ. Thus, to obtain the amplitude in units of δ and the frequency in units of F/mδ, we need only to specify ξ/δ and b/δ. The parameter b/δ is a free fitting parameter that will be set a value to get the best agreement with the numerical solution. In Fig. 3, we plot x 0 (t) obtained from the numerical solution of the time-dependent Gross-Pitaevskii equation together with the solution of Eq. (31) using the same initial conditions. The vortex is initially located at x 0 = −3.1δ from the Thomas-Fermi surface and then evolved starting from rest. The best agreement between the variational solution and the numerical one was obtained for b ≃ 5δ. The coherence length associated with this figure is calculated by fitting the numerical density profile to Eq. (8), using ξ as a fitting parameter. This was done for both the x-and y-cross sections of the vortex density profile as shown in Fig. 4. The average value of ξ turns out to be ξ ≃ 0.42δ. This is consistent with the Thomas-Fermi estimate ξ = 2mγ/h 2 ≃ 0.4δ, were we have used γ = gn(x 0 ) = −F x 0 . Since the coherence length obtained using the Thomas-Fermi approximation is approximately equal to the one obtained numerically, we conclude that during the whole time interval of Figs. 3 and 4, the vortex is moving in a region where the Thomas-Fermi approximation is valid. At first instance, we may take the average value of the coherence length, namely ξ = 0.4δ, and then substitute it in x eq and ω. This will lead to x eq ≃ −6δ. However, Fig. 3 indicates that the vortex oscillates around x eq ≃ −2.8δ. This discrepancy originates probably from the difference in geometry of the system in the variational calculation and the one in numerical calculation. In the variational calculation, the vortex is moving in an infinite background, but the energy integrations are performed within a cylinder of radius b = 5δ. Furthermore, the number of atoms within the cylinder is constant. In the numerical calculation, the system is a rectangular grid of dimensions 8δ × 40δ. The number of atoms within the whole grid is constant. Therefore, while in the variational calculation the number of atoms is kept constant within an area of πb 2 ≃ 80δ 2 , in the numerical calculation the number is kept constant within an area of 320δ 2 . Thus, there is a difference in the density which leads to a difference in the coherence length. In other words, the coherence length values shown in Fig. 4 are not exactly those which we should use in the variational calculation to seek agreement between the results of the two calculations. We found, as shown in Fig. 3, that the best agreement with the numerical values is obtained using ξ = 0.8δ.
We show also in Fig 3, vortex trajectories for initial distances x 0 (0) = −4δ and x 0 (0) = −5δ. There is a transient behavior in the beginning before acquiring the oscillatory behavior. For initial starting distances larger than 7δ, the oscillations disappear and the vortex moves in a straight path parallel to the surface as predicted by the variational calculation.

III. NUMERICAL SOLUTION
In this section, we describe our numerical procedure in solving the time-dependent Gross-Pitaevskii equation. Scaling length to δ, energy toh 2 /mδ 2 , and time to mδ 2 /h, Eq. (4) takes the following dimensionless form where the tilde denotes scaled quantities and |ψ| 2 is scaled to 1/4πaδ 2 . It is interesting to notice here that the scattering length disappears as a result of the scaling we employ here. We start by solving the time-independent Gross-Pitaevskii equation with a vortex-free initial state. We use the imaginary-time evolution method to solve this equation. The resulting densityñ 0 (x) is the background density profile for the vortex. We multiply this background density by a density profile of a vortex in a uniform background, namely χ 2 given by Eq. (8). The resulting density profile can not yet be used as an initial state for solving the time-dependent Gross-Pitaevskii equation. This is so since χ of Eq. (8) is not the exact density profile of the vortex. When we simulated the dynamics of the vortex using this density profile, the core of the vortex acquired breathing oscillations during the vortex motion. These oscillations produced circular density waves emitted out of the vortex core through the background and then reflected from the boundaries to interfere with vortex affecting its trajectory. Since we are interested only in the center-of-mass motion of the vortex, it is necessary to use the exact density profile of the vortex to avoid exciting such density waves. This can be obtained by evolving in imaginary time for a certain time period the densityñ 0 χ 2 . In other words, we resolve Eq. (37) with ñ 0 χ 2 as an initial wavefunction. The resulting profile contains the exact background and vortex core structure. This profile is then used as the initial profile for the time-dependent Gross-Pitaevskii equation. We found that density waves emerging from the vortex core have now almost disappeared. However, evolving in imaginary time results in the vortex sliding towards the surface of the condensate and ultimately disappearing there. In the present procedure, we evolve for only a finite time interval long enough to ensure adequate suppress of the density waves before the vortex reaches the surface. The distance from the surface reached by the vortex at the end of this time interval is to be considered as the initial distance for the real time evolution. We checked that our results for the real-time dynamics of the vortex are not affected by the length of this time interval once the density waves are significantly suppressed.
We employed a fourth-order Runge-Kutta method to integrate the time derivative of the time-dependent Gross-Pitaevskii equation. The spacial grids we used are 8δ × 40δ and 8δ × 80δ with spacial descretization of ∆x = 0.1δ and time descretization ∆t = 10 −5 . The long side of the grids is along the y-axis since the vortex is faster in that direction. We have used the boundary conditions that the derivative of density with respect to y vanishes at the two edges of the grid normal to the y-axis. The other two edges were left free. We have used larger grids to check that the boundaries do not affect the results obtained here and found that the above-mentioned grid sizes are adequate for that purpose. We have also calculated the total number of atoms and energy to find that they are fluctuating within less than 0.1%. The fluctuations were random and not correlated to the vortex motion.
To test the accuracy of our numerical calculation, we start by simulating the dynamics of a vortex near a hard wall in a condensate of uniform background density. The vortex velocity is calculated theoretically to be v y =h/2mx 0 [20]. In Fig. 2, we show that our code accurately generates such a dependence for the vortex velocity.
Having checked the accuracy of our numerical procedure, we simulate the dynamics of the vortex in the surface region of the condensate. In Fig. 5, we show the results in the form of a number of snapshots of the density of the condensate visualizing the realtime motion of the vortex core along the surface of the condensate. In Figs. 3 and 6, we plot the vortex coordinates x 0 and y 0 versus time. These two figures show that while the vortex is moving with a constant velocity in the y-direction, it oscillates in the x-direction as predicted in the previous section. To investigate the x 0 -dependance of v y , we calculate the numerical value of v y for a number of different values of x 0 . The results are shown in Fig. 2 where we also plot the theoretical prediction of Anglin [20] and the result of our variational calculation. The formula derived by Anglin is supposed to be valid away from the Thomas-Fermi surface. Figure 2 shows indeed that while this formula agrees well with the numerical calculation for large x 0 , it deviates from it for values of x 0 of order δ. On the other hand, this figure shows that our formula, derived in subsection II B, is in a better agreement with the numerical data especially for x 0 < 10δ where the formula of Anglin starts to deviate from the numerical values.
We have also calculated the coherence length along both the x and y-directions, which is of the order of the size of the core of the vortex. This is plotted in Fig. 4. This figure shows that the core sizes change slightly over the whole time interval considered. The average value of ξ in this figure agrees with the Thomas-Fermi estimate as explained in subsection II C. In the previous section, we performed the calculation using the Thomas-Fermi approximation for the background density and neglected the dynamics of the size of the core. It is clear now that this is indeed justified here.

IV. CONCLUSIONS
We have shown that our variational calculation describes at least qualitatively the motion of a vortex near the surface of a Bose-Einstein condensate, and generates results that agree with previous calculations [16,19]. These qualitative results were then supported and quantified using a numerical simulation of the vortex motion.
One of the main findings of this work is Eq. (20) which gives the y-component (parallel to the surface) of the vortex velocity as a function of it is distance from the Thomas-Fermi surface x 0 . It agrees qualitatively with that of Anglin [20], and describes more accurately v y for values of x 0 of order of δ as shown in Fig. 2. This figure shows clearly that Eq. (20) can be used for both small and large x 0 /δ, unlike the result of Anglin [20] which is accurate only for large x 0 /δ.
Another main result of this work is that, if the vortex is moving within few coherence lengths ξ from the Thomas-Fermi surface, it will have an oscillatory component of its velocity in the direction normal to the surface. These oscillations disappear for x 0 ≪ −δ and the vortex will have only the y-component of the velocity. To further understand the physics of these oscillations, we have shown in subsection II B that the vortex feels an effective potential that is a function of x 0 and has local minimum for negative values of x 0 . This effective potential is composed mainly of a sum of a mean-field-energy and a trapping potential energy part as shown in Eq. (22). These two parts compete such that the effective potential has a minimum at x eq = −(δ/ξ) 2 δ. Using a collective coordinates approach, we have shown that the vortex oscillates around this minimum with a frequency ω = (F 4 ξ 2 /8mγ 3 ) ln (b/ξ). As discussed in subsection II B, this frequency is of order less than |x 0 | −4 ln (|x 0 |/ξ). Both this frequency and the x 0 -dependent part of U (x 0 ) vanish for large x 0 . This means that the vortex oscillations will disappear away from the surface of the condensate and will be noticeable only when ξ/x 0 is not negligible. To the best of our knowledge, these oscillations are predicted here for the first time. As a possible alternative explanation to these vortex oscillations near the surface, we mention the possibility that the vortex being dragged by a collective or surface wave. In our numerical procedure, the evolution in imaginary time which generates the initial state for the real time evolution, may include some excitations of collective modes. Further investigation is needed to verify this argument.
It was shown in Ref. [20] that the energy of a vortex near the Thomas-Fermi surface is not infrared divergent as in the bulk. This is so since the vortex is localized to within its distance from the surface. To avoid these infrared divergencies in the energy, we have introduced an infrared cut-off on the energy. This is equivalent to using the infrared-free flow pattern of Ref. [20] which gives finite energy.

APPENDIX A: DETAILS OF THE CALCULATION OF THE LAGRANGIAN
In this Appendix we present the details of calculating the lagrangian Eq. (16). This is performed by substituting the trial trial wave function Eq. (9) in Eq. (11).
We start by calculating the energy functional Eq. (13). According to the geometry illustrated in Fig. 1, the cartesian coordinates x ′ = x−x 0 and y ′ = y −y 0 with origin at the core of the vortex are related to the polar coordinates ρ(x ′ , y ′ ) and φ(x ′ , y ′ ) by ρ(x ′ , y ′ ) = x ′ 2 + y ′ 2 and φ(x ′ , y ′ ) = tan −1 (y ′ /x ′ ). This leads to that ∂φ/∂x = ∂φ/∂x ′ = cos φ/ρ and ∂φ/∂y = ∂φ/∂y ′ = − sin φ/ρ. In addition, we have ∂ρ/∂x = ∂ρ/∂x ′ = sin φ and ∂ρ/∂y = ∂ρ/∂y ′ = cos φ. The derivative of the trial wave function ψ(x, y, t) with respect to x, can thus be shown to take the form The contribution to the kinetic energy is obtained by integrating |∂ψ/∂x| 2 with respect to ρ and φ. The latter integral reduces to 2π 0 dφ ∂ψ ∂x where χ ρ denotes dχ(ρ)/dρ. It should be noted that we neglected here the contributions of order (ρ/x 0 ) 2 (ln (αρ/x 0 )) 2 that originate from the logarithmic term of the phase Φ(ρ, φ). In deriving this term, Anglin [20] performed a perturbative expansion of the phase of the vortex near the core using ξ/x 0 as the perturbation parameter, where ξ is the coherence length defined in subsection II A. Since in Anglin's calculation only terms of order (ξ/ρ) ln (αρ/x 0 ) are kept, it will be inconsistent if we keep here higher order terms. We keep, however, the higher order terms that do not originate from the logarithmic term of the phase. Similarly, for the derivative with respect to y, we have The kinetic energy is then given by The potential energy associated with the external force is The mean-field energy is Combining these three contributions, we find the energy functional per atom where γ = gN/N 1 l and we have defined the coefficients The relevant energy is the energy change due to the presence of the vortex. Therefore we have to subtract from the above energy functional the energy of a vortex-free background Thus, the energy difference associated with vortex reads The time derivative of the trial wave function is given by Using that ∂ρ/∂t = (∂ρ/∂x ′ )(∂x ′ /∂t) + (∂ρ/∂y ′ )(∂y ′ /∂t) = −ẋ ′ 0 sin φ −ẏ ′ 0 cos φ and similarly ∂φ/∂t = (−ẋ ′ 0 cos φ + y ′ 0 sin φ)/ρ and Eq. (7), the part of the lagrangian (per atom) containing the time derivatives takes the form The lagrangian per atom, Eq. (16), is then obtained.

ACKNOWLEDGMENTS
I would like to thank Halvor Nilson from NORDITA for useful discussions. I would also like to thank the anonymous referee whose comments have significantly enhanced this paper. The planner geometry and coordinate system used in this paper. The small circle represents the core of the vortex which is of the order of the coherence length ξ. The coordinates of the core of the vortex are x0 and y0. The radius of the condensate is R. The bulk of the condensate is located in the region x < 0. The point P has cartesian coordinates x ′ and y ′ related to the polar coordinates ρ and φ by x ′ = ρ sin φ and y ′ = ρ cos φ. The vortex velocity in the y-direction as a function of the magnitude of the distance between the vortex core and the Thomas-Fermi surface. The velocity is in units of δ/τ . The surface depth δ is defined byh 2 /mδ 2 = F δ, and the time τ is defined by τ =h/F δ. The dots on the upper curve are the results of the numerical solution of a vortex moving near the surface of the condensate. The dots on the lower curve correspond to a vortex moving parallel to a hard wall in a uniform density background. The solid curve is the result of the variational calculation of section II, given by Eq. (20). The long-dashed curve corresponds to the expression derived by Anglin [20]. The short-dashed curve is given by vy =h/2mx0.  The right column is a side view in which the y-axis is normal to the picture towards the viewer. The vertical line in the right column of pictures is to reference the core location with respect to the initial position, and shows clearly that the vortex core oscillates in the x-direction.