Surface excitations, shape deformation and the long-time behavior in a stirred Bose-Einstein condensate

The surface excitations, shape deformation and the formation of persistent current for a Gaussian obstacle potential rotating in an highly oblate Bose-Einstein condensate(BEC)are investigated. Vortex dipole can be produced and trapped in the center of the stirrer even for slow motion of the stirring beam. When the barrier angular velocity is above some critical value, the condensate shape can be deformed remarkably according to the rotation frequency due to the existence of plenty of surface wave excitations. After a long enough time, a few vortices are found to be left either trapped in the condensate or pinned by the obstacle, a vortex dipole or several vortices can be trapped at the beam center, which enables the possibility of vortex manipulation.


Introduction
Quantized vortex is a topological singularity in a superfluid or superconductor, where the phase of the order parameter varies by an integer multiplying 2π when following a closed path around the defect. Due to its topological nature and the conservation of circulation, a vortex can be eliminated only by annihilation with an antivortex or moving to the boundary of the system [1].
One of the important developments in atomic BEC is the recent experimental observations of a vortex dipole when a repulsive Gaussian obstacle moves through a pancake-shaped condensate [38,39]. Subsequently, much theoretical and experimental effort has been directed toward the topic of vortex dipoles. The studies mainly focus on the vortex shedding mechanism [40][41][42][43], topological excitation [42,44], and persistent current [45]. Additionally, vortex dipoles induced in oscillating potential [46], at finite temperatures [47], and under spin-orbit coupling [48] have been extensively studied. The experimental technique using two blue-detuned laser beams as "tweezers" provides realistic possibilities to pin and manipulate vortices [49][50][51]. Recent progress in a toroidal geometry [52][53][54][55][56][57][58] has opened a new prospect to study the superfluidity and vortex excitation using a rotating barrier or weak link [53,59,60]. When a weak link rotates in an annular BEC, it is found that surface modes can be excited, and the vortex can come into the condensate from the outside [57]. Furthermore, a persistent flow in a toroidal trap can be created by stirring with a rotating barrier [54]. For a highly oblate BEC stirred by a laser beam, studies have been done on the shape deformation elliptically and triangularly [61,62]. However, there is a lack of extensive study on this topic, and, more importantly, there have been no investigations on shape deformation for larger angular momentum, such as when the quantum number is 6 and 7. Moreover, the vortex dipoles inside the obstacle on their nucleation, the splitting properties, and the long-time dynamical behavior is less concerned and studied but is of great interest.
In this paper, by solving the time-dependent damped Gross-Pitaevskii (GP) equations, we reveal the shape deformation and long-time dynamical behaviors of a highly oblate BEC stirred by a rotating laser beam. The shape of the condensate deformed heavily due to excitations of many surface waves with relatively lower angular momenta. After a long enough time, a small numbers of vortices can be left either trapped in the condensate or pinned by the obstacle. On the other hand, stirring BEC by a rotating blue-detuned laser beam can cause a vortex dipole or several vortices to be pined at the beam center.

Model
Consider a single-component BEC described by the normalized macroscopic wave function ψ(r, t).
In the mean-field framework, the dynamics of a system with N weakly interacting identical atoms close to thermodynamic equilibrium and being subjected to weak dissipation can be described by the damped GP equation [63]: where V(r) = 1 2 m(ω 2 t x 2 + ω 2 t y 2 + ω 2 z z 2 ) is the axially symmetric harmonic trap potential, and ω t , ω z are the radial and axial trap frequencies. The parameter γ in the GP equation is a dimensionless, phenomenological, damping constant. It takes into account the quantum and thermal fluctuations from the background and is introduced to fit the experiment. In studies of vortex excitations in an annular BEC [56], γ is chosen to be 0.0015, by which the measured lifetime of the vortices is consistent with the experimental result in [54]. The damped GP equation has been extensively employed to study the dynamics of systems in the presence of thermally induced dissipation [24,50,51,56,57].
We focus in this paper on a highly oblate BEC with ω z ω t . In this extreme limit, the axial dimension is tightly suppressed that the motion along the z direction can be neglected and atoms can move only within the x − y plane. The normalized ψ can thus be written as ψ(r, t) = ψ 2D (x, y, t)φ 0 (z), where φ 0 (z) = (πa 2 z ) −1/4 e −z 2 /2a 2 z is the ground state for the potential V z (z), with a z = √h /mω z the characteristic length of the vertical harmonic oscillator. The atom-atom contact interaction is g = 4πh 2 a s /m, with a s as the s-wave scattering length. In the numerical computations, we discretize the x − y plane into a square lattice with x(i) = ia and y(j) = ja. The artificial lattice constant a must be much less than the characteristic length l = √h /mω t of the axial harmonic oscillator to validate the simulations [64]. Furthermore, a should also satisfy the condition of a ξ, with ξ as the condensate healing length. For the time evolution, the fourth-order Runge-Kutta method is employed at each time step. The central-difference formula is used mainly to calculate the kinetic term. Introducing dimensionless ψ(i, j), by substituting ψ(r) with 1 √ a 2 a z ψ ( i, j), we thus obtain the following lattice-version GP equations: where t 0 =¯h 2 2ma 2 ,Ṽ = 1 2 mω 2 t a 2 , andg = 4πa sh 2 ma 2 a z . Since all the parameters t 0 ,Ṽ,g, andh/t have the scale of energy, it is convenient to introduce dimensionless parameters V = V/t 0 , g = g/t 0 , and t = t/(h/t 0 ), measured in units of t 0 . All of these parameters are actually only dependent of a/l and a s /a z . A straightforward analysis leads to the following expressions of V = ( a l ) 4 and g = 8π a s a z . Note that g is essentially independent of the artificial lattice constant a.
We assume the system consists of 87 Rb atoms. Thus, we have m ≈ 87m p , with m p as the mass of proton. We choose the trap frequency ω t = 2π × 10 Hz. l is then estimated to be about 3.4 µm.
When the square lattice we study takes a typical size of 200 × 200 and (a/l) 2 is chosen to be 0.008, the system has a size of about 60 µm × 60 µm andhω t /t 0 = 2(a/l) 2 = 0.016,h/t 0 = 0.2 ms, and V = 0.000064. In addition, to guarantee both the convergence and efficiency of iteration of the GP equations, dt is chosen to be between 10 −4 and 10 −2 in numerical calculations. We control the convergence from the density distribution until the precision reaches the level of 10 −10 in each mesh.

Results and Discussions
When a blue-detuned laser beam is rotating uniformly with angular velocity Ω in BEC, the obstacle produced by it can be well described as a moving Gaussian potential: where {x 0 (t), y 0 (t)} = {R cos(Ωt), R sin(Ωt)}, with R as the distance between the obstacle and trap center, and σ and V 0 are the laser beam waist and barrier height of the potential, respectively. The detailed behavior of the stirred BEC depends sensitively on σ/ξ and V 0 [41][42][43]49,53,65], and in the following they are fixed to be V 0 = 8t 0 , σ = 3.5ξ. In such a quasi-2D atomic system, at the trap center, the chemical potential µ is estimated to be µ = n 0 g ≈ 0.5t 0 , where n 0 is the atomic density |ψ(i, j)| 2 .
Therefore, the condensate healing length ξ =¯h √ 2mµ and the sound velocity v s = µ m are on the order of 0.2l and 1.15 mm/s, respectively. By starting from the ground state of the BEC in the presence of a static obstacle, the uncontrollable excitations are prevented from an abrupt motion of the obstacle. Before the obstacle moves with a constant angular velocity Ω, the angular velocity is assumed to increase linearly with time from zero until reaching Ω. In the following, γ is set to be 1.5 × 10 −3 . We find our results do not depend qualitatively on the specific value of γ when it is in the range of 0.001-0.005.
According to the Landau criterion, a superfluid becomes dissipative when it flows above a certain critical velocity v c via generating its elementary excitations such as phonons and vortices. For a homogeneous system, the Landau critical velocity is equal to the speed of sound v s [66]. The critical velocity is believed to provide an upper bound for a moving hard cylinder in a uniform BEC and depends on the barrier height V 0 and width σ in a harmonic tap [41]. Here, the critical velocity v c is defined as the value of the obstacle velocity v = RΩ, at which the vortex excitation begins to be generated into the BEC by the obstacle and depends on the distance R from the BEC center. However, even when v is much smaller than v c , a vortex dipole can be produced and tightly trapped by the obstacle. Since the obstacle potential is finite in spite of being large, the local density n of the BEC within the small space occupied by the obstacle is rather small but finite. The small local density thus lowers the speed of sound locally, leading to the local excitations within the obstacle according to the Landau criterion. When the obstacle velocity v is below v c but not too small, the vortex dipole direction can be identified from the phase profile of the BEC, which rotates uniformly within the obstacle but always remains perpendicular to the motion direction of the obstacle, as exhibited in Figure 1. Besides the critical velocity v c , there exists a critical angular velocity Ω c = min{ω l /l} according to an analog of the Landau criterion [67][68][69], where ω l is the surface-wave excitation energy. For smaller l, by solving the hydrodynamic equations of a superfluid, one can obtain the general dispersion law ω l = √ lω t [70]. For the parameters chosen, Ω c is estimated to be about 0.3ω t , above which surface waves can be excited. However, the surface excitations with larger l can hardly be identified from the density or phase profiles of the BEC obtained numerically. On the other hand, for smaller l, ω l /l = ω t / √ l can also be viewed as the minimum excitation frequency for generating a surface excitation with angular momentum l, indicating that the surface excitations with smaller l can be possibly found upon increasing Ω up to a sufficient large value. It is found in our calculations that, when Ω = 0.6ω t , the density of the BEC is deformed triangularly. Figure 1(a2)-(a3) also show the similar behavior of the condensate when the rotating frequency Ω is slightly larger than ω t /2(ω t / √ 5), where the condensate is deformed tetragonally(pentagonally). Furthermore, by investigating the time evolution of density, we find that the shape deformation of the condensate is rotating rigidly nearly at the same angular velocity as the obstacle, as can be seen in Figure 1, where the two density profiles in each column show the evolution of the density at two instants. This feature can be explained as follows. The condensate can be viewed to be composed of two parts: the surface and core parts. The surface part with strong shape deformation consists of atoms which carry the surface excitations and their angular momenta, while the core part is left nearly motionless. The surface excitations with angular momentum l produce a density deviation δn ∝ e i(lθ−ω l t) = e il(θ− ω l l t) , indicating that the corresponding shape deformation is rotating with an angular velocity ω l /l, while the rotating obstacle can be seen as a density perturbation δn ∝ δ(θ − Ωt) = ∑ 1 2π e in(θ−Ωt) . The resonance occurs when Ω approximately equals ω l /l, where many surface waves with angular momentum l are excited, leading to a heavy shape deformation in the form of an l-regular polygon. In this situation, the surface part is also rotating nearly synchronously with the obstacle.
On the other hand, the surface part is believed to contain much fewer atoms than the core part, so the average angular momentum per atom is approaching a relative small value (about 0.15 for the case in Figure 1(c2) in the long-time evolution). Notice that the contribution to angular momentum from the vortex dipole is expected to be quite small. As shown schematically in Figure 1d, the vortex dipole's contribution to the angular momentum is estimated to be ∼ d 0 R TF 0.1 with d 0 the vortex-antivortex distance and R TF the TF radius, since d 0 is found to be much less than the healing length ξ.
Since Rω t / √ 2 > v c for R = 1.8l, the elliptical deformation due to l = 2 surface waves cannot be observed. However, for a smaller R = 1.3l, our calculations confirm the existence of the elliptically deformed BEC for Ω 0 = ω t / √ 2, which is consistent with the experimental observations [5,67]. In Figure 2, we plot our results of the critical angular velocity for surface wave excitations when l varies from 2 to 8 and comparing them with the analytical value of ω t / √ l.  We now turn our attention to the critical velocity and mechanism of the vortex nucleation. When a circular cylinder is passing through a two-dimensional homogenous system, the critical velocity is calculated numerically and estimated to be v c = 0.37v s in the large-cylinder limit R ξ, with R as the radius of the cylinder [71]. Experimentally, the minimal measured value of v c /v s for vortex shedding in a highly oblate BEC reaches 0.1 [41]. In our model, when R = 1.8l, v c is estimated to be v c = 0.21v s | r=R = 242 µm/s, at which the vortex dipole starts to split. This critical value can also be confirmed from the drag force defined by F = ih∂ t dr(Ψ † ∇Ψ) [17]. The drag forces along the x and y directions at 0.20v s and 0.21v s are illustrated in Figure 3. Under the same velocity, both x and y directions gradually exhibit periodic oscillations with the same frequency due to the periodic motion of obstacle. While at v c the drag forces gradually decrease since the accumulated energy is released [72]. The antivortex starts to separate from the vortex and move toward the barrier edge, while the vortex is still pinned at the obstacle center, as schematically shown in Figure 4a. Meanwhile, a sharp corner of the density hole appears (Figure 4b), long before the antivortex leaves the barrier edge. The antivortex is then released from the obstacle and moves gradually to the boundary of the condensate. When v is slightly above 0.26v s , after the first antivortex escapes from the obstacle, another vortex-antivortex pair will be generated and subsequently the new antivortex starts to escape, leaving the two vortices trapped at the center by the obstacle (see Figure 4). When v > 0.3v s , more vortex-antivortex pairs will be generated during the whole stirring process and more vortices can be trapped by the obstacle. It is interesting to compare the microscopic mechanism of vortex nucleation to the case of a weak link rotating in a torodial BEC. As was shown in [57], at a low rotation rate above the threshold value, a vortex enters the weak link from the outside and approaches the antivortex to create a vortex-antivortex pair, since the critical velocity for vortex nucleation is reached first at the outer edge of the condensate.  After stirring for over 100 circles, the BEC reaches a metastable state, which shows nearly time-independent features. This means the obstacle will stop shedding vortices after a long-time stirring. This can be better understood in the limit of V 0 being infinite. The vortices trapped by the obstacle or the condensate have changed the distribution of the superfluid velocity around the obstacle. According to Landau criterion, no vortex will be excited, if all velocities are below the critical value, leading to the formation of persistent current after a long-time stirring. Now we study the long-time behavior of the BEC and focus on the vortex number left in the condensate. After a long enough time, all antivortices will leave the condensate, and only vortices are left. Among the remaining vortices, several are pinned by the obstacle, and the others are loosely trapped near the trap center. The nearly time-independent angular momentum also confirms the stability of the long-time behavior of the BEC. In the parameter region we studied, three vortices at most can be trapped by the obstacle. These results are demonstrated in Figure 5. The number of trapped vortices is summarized in Table 1. If the laser beam stops moving suddenly, the trapped vortices will be released, but finally one or two of them will be trapped at the obstacle center and stabilized to form persistent current state. When the laser beam is ramping off, the trapped vortices will become free and will then move to the boundary. Comparing with the method for vortex manipulation using the chopstick beams [50], the rotating obstacle in our model plays a key role in the mechanism of pinning vortices, which provides another way to manipulate a vortex. For the parameters we studied here, we only focus on the velocity regime v < 1.7v c , since, when v > 2v c , many vortices are generated quickly such that the stirred BEC is in a turbulent state. Within the regime v c < v < 1.7v c , it is not from the boundary but within the obstacle center in the form of a vortex pair that the vortices are created.   Figure 5. The long-time dynamical behavior of the condensate when the obstacle velocity is above the critical value v c with v c = 255 µm/s at R = 3.6l and v c = 242 µm/s at R = 1.8l, by stirring the BEC for over 100 circles. The upper (lower) two panels denote the density and phase profiles with R = 3.6l (R = 1.8l). The middle panels denote the angular momenta per atom corresponding to the cases from (a1) to (a4), respectively.

Conclusions
We have performed numerical calculations of the quasi-two-dimensional GP equation to investigate the surface wave excitations and the long-time behavior of vortices in a stirred, highly oblate condensate. Surface waves can be excited with the condensate shape being deformed heavily at the corresponding rotation frequency. A vortex-antivortex pair can be created in the obstacle center even when the obstacle velocity is relatively small. Once the obstacle velocity reaches a critical value, the antivortex starts to separate from the vortex and then leaves the obstacle regime. Furthermore, after a long enough time, a small number of vortices are found to be left either trapped in the condensate or pinned by the obstacle, and the number of them depends on the velocity and position of the obstacle.

Conflicts of Interest:
The authors declare no conflict of interest.