Vortex Rings and Lieb Modes in a Cylindrical Bose-Einstein Condensate
S. Komineas1 and N. Papanicolaou2
Physikalisches Institut, Universit¨ at Bayreuth, D-95440 Bayreuth, Germany Department of Physics, University of Crete, and Research Center of Crete, Heraklion, Greece emails: email@example.com firstname.lastname@example.org (Dated: February 1, 2008)
We present a calculation of a solitary wave propagating along a cylindrical Bose-Einstein trap, which is found to be a hybrid of a one-dimensional (1D) soliton and a three-dimensional (3D) vortex ring. The calculated energy-momentum dispersion exhibits characteristics similar to those of a mode proposed sometime ago by Lieb within a 1D model, as well as some rotonlike features.
arXiv:cond-mat/0202182v1 11 Feb 2002
PACS numbers: 05.45.Yv, 03.75.Fi, 05.30.Jp
Approximately forty years ago, Lieb and Liniger  employed the Bethe Ansatz to obtain an exact diagonalization of the model Hamiltonian that describes a onedimensional (1D) Bose gas interacting via a repulsive δ function potential. Based on this solution Lieb  proposed an intriguing dual interpretation of the spectrum of elementary excitations, either in terms of the familiar Bogoliubov mode, or in terms of a certain kind of a particle-hole excitation which will hereafter be referred to as the Lieb mode. The energy-momentum dispersions of the Bogoliubov and the Lieb modes coincide at low momenta, where they both display the linear dependence characteristic of sound-wave propagation in an interacting Bose-Einstein Condensate (BEC), but the Lieb dispersion is signi?cantly di?erent at ?nite momenta where it exhibits a surprising 2π periodicity. Yet this new mode remained somewhat of a theoretical curiosity because of the absence of a physical realization of a strictly 1D Bose gas. In a curious turn of events, the above subject reemerged in connection with the solitary waves calculated analytically within a 1D classical Gross-Pitaevskii (GP) model [3, 4], which were later shown to provide an approximate but fairly accurate description of the quantum Lieb mode for practically all coupling strengths [5, 6]. More importantly, the same solitary waves motivated experimental observation of similar coherent structures in BEC of alkali-metal atoms through phase imprinting or phase engineering [7, 8]. Therefore, an opportunity presents itself for experimental realization of the Lieb mode. It is clear that a theoretical investigation based on effective 1D models often used to describe cylindrical traps will lead to the prediction of a Lieb mode, by arguments similar to those employed in the original 1D classical model [5, 6]. However, one should also question the stability of the corresponding solitons within the proper 3D environment of a realistic trap [9, 10]. For example, dark solitons created within a spherical trap were observed to decay into vortex rings . This observation was also supported by a numerical calculation in which an initial dark-soliton con?guration is let to evolve in time according to the GP equation. Therefore, although the production of solitary waves
through phase imprinting clearly suggests some distinct 1D characteristics, it should be expected that a proper understanding of such waves must also account for 3D e?ects that are inevitably present in a realistic BEC. One could envisage a picture in which the actual solitary wave is a hybrid of a 1D soliton and a 3D vortex ring. It is the aim of the present paper to make the above claim precise by calculating solitons propagating along a cylindrical trap without making a priori assumptions about their e?ective dimensionality. Our approach was motivated by the calculation of vortex rings in a homogeneous BEC due to Jones and Roberts  and a similar calculation of semitopological solitons in planar ferromagnets . Explicit results will be presented for a set of parameters that roughly correspond to the experiment of Ref. . Thus we consider a cigar-shaped trap ?lled with 87 Rb atoms, with a transverse con?nement frequency ω⊥ = 2π × 425 Hz and a corresponding oscillator length a⊥ = ? h/mω⊥ ≈ 0.5?m. The coupling constant is written as Uo = 4π ? h2 a/m where a ≈ 50 ? A is the scattering length. We make the approximation of an in?nitely-long cylindrical trap with average linear density ν = 0.2 atoms/? A and introduce the dimensionless combinations of parameters γ = νa = 10 and γ⊥ = νa⊥ = 103 . Finally, rationalized units are de?ned √ through the rescalings t → t/ω⊥ , r → a⊥ r, and Ψ → ν Ψ/a⊥ . The energy functional extended to include a chemical potential is then given by W = 1 2 ?Ψ? ?Ψ + ρ2 Ψ? Ψ + g (Ψ? Ψ)2 ? 2? Ψ? Ψ dV ,
(1) where g = 4πγ and ρ2 = x2 + y 2 , and yields energy in units of γ⊥ ? hω⊥ while the chemical potential ? is still measured in units of h ? ω⊥ . Similarly, the conserved linear momentum along the z axis given by the usual de?nition P= 1 2i Ψ? ? Ψ ? Ψ? ? Ψ ?z ?z dV = n ?φ dV , ?z (2)
is measured in units of h ? ν = γ⊥ (? h/a⊥ ). In the second step of Eq. (2) we employed the usual hydrodynamic variables √ de?ned from Ψ = n eiφ . An important ?rst step in out calculation was to obtain accurate information about the ground state and the cor-
FIG. 1: Radial dependence of the solitary wave function at v = c/2 for various positive values of ξ in steps of δξ = 0.1. The result for negative ξ is obtained through the symmetry relations Re Ψ(ρ, ξ ) = Re Ψ(ρ, ?ξ ) and Im Ψ(ρ, ξ ) = ?Im Ψ(ρ, ?ξ ). Distances are measured in units of a⊥ .
FIG. 2: Radial dependence of the local particle density n = |Ψ|2 for the four special cases discussed in the text. Conventions are the same as in Fig. 1.
responding linear (Bogoliubov) modes. The ground-state wave function Ψ = Ψ0 (ρ) was calculated as a stationary point of the energy functional (1) by a variant of a relaxation algorithm  and is normalized according to ∞ 2 0 2πρ dρ |Ψ0 | = 1 to conform with our choice of rationalized units. The chemical potential was found to be ? = 6.4324 for γ = 10. We have also recalculated the lowest branch of the Bogoliubov spectrum from which we extracted the speed of sound c = 1.77 in units of a⊥ ω⊥ . This numerical estimate at γ = 10 is consistent to within 1% with the Thomas-Fermi approximation c = γ 1/4 which was previously derived in a number of papers [15, 16, 17] and was employed for the analysis of experimental data . Thus we turn to the calculation of axially-symmetric solitary waves described by a wave function of the form Ψ = Ψ(ρ, ξ ) where ξ = z ? vt and v is the constant velocity along the z axis. Such a wave function satis?es the stationary di?erential equation ? iv ?Ψ 1 1 = ? ?Ψ + ρ2 Ψ + g (Ψ? Ψ)Ψ ? ?Ψ , (3) ?ξ 2 2 ? = 1 ? ?2 ?2 + , + ?ρ2 ρ ?ρ ?ξ 2
solution of Eq. (3) must satisfy the virial relation 1 ? Ψ? ? Ψ g + ρ2 Ψ? Ψ + (Ψ? Ψ)2 ? ?Ψ? Ψ dV , 2 ?ξ ?ξ 2 (4) obtained by a standard scaling argument. Solutions of Eq. (3) were obtained by a sophisticated Newton-Raphson iterative algorithm [12, 13] which will not be described here in any detail. For values of the velocity v near the speed of sound c, the calculated solitary wave is a rarefaction pulse that may be thought of as a mild soundlike disturbance of the ground state. The dominant features of the solitary wave are pronounced as the velocity is decreased to lower values and become reasonably apparent at v = c/2, as shown in Fig. 1 which provides a complete illustration of the wave function through its real and imaginary parts. A partial but more transparent illustration is given in frame I of Fig. 2 where we depict the radial dependence of the local particle density n = |Ψ|2 for various values of ξ . It is clear that the density near the center of the soliton (ξ = 0) vanishes on a ring with a relatively large radius R ? 2.8, thus a vortex ring is beginning to emerge. The features of the vortex ring become completely apparent, and its radius is tightened, as we proceed to yet smaller values of the velocity. A notable special case is the static (v = 0) vortex ring with radius R ? 1.8 illustrated in frame II of Fig. 2, which is far from being a completely dark (black) solitary wave. One would think that pushing the velocity v to negative values would retrace the calculated sequence of solitary waves backwards. In fact, our algorithm continues vP =
which is supplemented by the boundary conditions that Ψ vanish for ρ → ∞, and |Ψ| → |Ψ0 (ρ)| for ξ → ±∞. The phase of the wave function is not ?xed a priori at spatial in?nity, except for a mild restriction implied by the Neumann boundary condition ? Ψ/?ξ = 0 imposed at ξ → ±∞ mainly for numerical purposes. Finally, a
3 to converge to vortex rings of smaller radii until a critical velocity v = ?v0 ? ?0.475 c is encountered where the vortex ring achieves its minimum radius Rmin ? 0.8 and ceases to exist for smaller values of v . The terminal state at v = ?v0 is illustrated in frame III of Fig. 2. We have thus described a sequence of solitary waves with velocities in the range ?v0 < v < c, which does not contain a black soliton. An equivalent sequence is obtained in the range ?c < v < v0 simply by reversing the relative sign between the real and imaginary part of the wave function, as is evident from Eq. (3). It is now important to calculate the energy-momentum dispersion. The excitation energy is de?ned as E = W ? W0 where both W and W0 are calculated from Eq. (1) applied for the solitary wave Ψ(ρ, ξ ) and the ground state Ψ0 (ρ), respectively. The presence of the chemical potential in Eq. (1) provides the compensation that is necessary in order to compare energies of states with the same number of particles. Similarly, the relevant “physical” momentum is not the linear momentum P of Eq. (2) but the “impulse” Q de?ned in a manner analogous to the case of a homogeneous gas [6, 12]: Q = δφ ≡ (n ? n0 )
3 IV γ = 10
II 2 III E I 1
?φ dV = P ? δφ ?z
FIG. 3: Energy E in units of γ⊥ (? hω⊥ ) versus impulse Q in units of γ⊥ (? h/a⊥ ). The solid line corresponds to the main sequence of solitary waves discussed in the text, and the dashed line to the auxiliary sequence that contains a black soliton (point IV).
(5) Fig. 3 was calculated on the basis of the sequence of solitary waves with velocities in the range ?v0 < v < c. The branch for ?c < v < v0 may be thought to correspond to negative Q, as usual. However, there is an implicit 2π periodicity because of the appearance of the angle φ in the de?nition of the impulse in Eq. (5). Therefore, a more natural representation of the spectrum seems to be obtained by restricting Q to the “?rst Brillouin zone” and thus completing Fig. 3 by its mirror image around Q = π. We have also carried out a detailed numerical calculation of the Bogoliubov dispersion for comparison. But the main point can be made here by assuming the model dispersion ω = |q |(c2 + q 2 /4)1/2 where the frequency ω is measured in units of ω⊥ and the wave number q in units of 1/a⊥. Translated into the units employed in Fig. 3, the model Bogoliubov dispersion reads 2 2 2 ≈ 106 E = |Q|(c2 + γ⊥ Q /4)1/2 where a huge factor γ⊥ makes its appearance. Therefore, although the Bogoliubov and Lieb dispersions coincide at low Q, they diverge quickly at the scale of Fig. 3, a notable di?erence from the homogeneous 1D model [2, 6]. In other words, the two modes operate at rather di?erent energy and momentum scales in a cylindrical trap. The last question we address is whether or not there exists an independent sequence of solitary waves that reduces to a black soliton at v = 0. We used a model blacksoliton con?guration as input in our iterative algorithm which converged to the true static soliton illustrated in frame IV of Fig. 2 that is indeed a black soliton. We then incremented the velocity to either positive or negative values with symmetrical results. For de?niteness, we follow the solution to negative v and ?nd that it ex-
2πρ dρ n0 (ρ)[φ(ρ, z = ∞) ? φ(ρ, z = ?∞)] ,
where n0 = |Ψ0 (ρ)|2 is the ground-state particle density and δφ is the weighted average of the phase di?erence between the two ends of the trap. Here we simply postulate the validity of de?nition (5) and note that the corresponding group-velocity relation v = dE/dQ was satis?ed to an excellent accuracy in our numerical calculation. On the other hand, the virial relation (4) was veri?ed using the standard linear momentum P of Eq. (2), as expected. The calculated dispersion E = E (Q) is depicted in Fig. 3 by a solid line along which we place the symbols I, II and III that correspond to the special cases of the solitary wave discussed above. At low Q the dispersion is linear, E = c |Q|, a feature that it shares with the Bogoliubov dispersion. However, the energy now reaches a maximum at point II where v = 0 and Q is slightly greater that π . Interestingly, the group velocity becomes negative in the region (II,III) or, equivalently, the impulse is opposite to the group velocity. This rotonlike behavior could develop into a full-scale roton if the terminal point III actually turns out to be an in?ection point beyond which the group velocity may begin to rise again. We have not been able to somehow continue our sequence of solitary waves beyond point III, nor have we ruled out such a possibility. In any case, the picture just described, together with the fact that the calculated radius of the vortex ring is monotonically decreasing along the sequence (I,II,III) comes close to the Onsager-Feynman view of a roton as the ghost of a vanished vortex ring . The branch of the dispersion shown by a solid line in
4 ists only until one encounters the same critical velocity ?v0 discussed earlier in the text. The corresponding terminal state is precisely the same with the one reached through our original sequence of solitary waves and illustrated in frame III of Fig. 2. For intermediate values of v the solution is again a vortex ring with constant radius at R = Rmin ≈ 0.8. Accordingly, the calculated energymomentum dispersion shown by a dashed line in Fig. 3 joins the original sequence through a cusp at the terminal point III, a situation that is reminiscent of the calculation of Jones and Roberts . The same authors together with Putterman later argued  that the solitary waves that correspond to the upper branch are actually unstable, a conclusion that might be valid in the present case as well. A summary view of our results is given in Fig. 4. We have thus constructed a family of solitary waves which exhibit some quasi-1D features, such as the appearance of a nontrivial phase di?erence δφ that is important for phase engineering [7, 8], but are otherwise bona?de 3D vortex rings.
FIG. 4: Contour levels of the local particle density n in a plane that contains the z axis and cuts across the cylindrical trap. The complete 3D picture may be envisaged by simple revolution around the z axis. Regions with high particle density are bright while regions with zero density are black. The four special cases considered are the same as in Fig. 2.
We are grateful to G.M. Kavoulakis for numerous discussions that guided us through the extensive recent literature on BEC. S.K. acknowledges ?nancial support from the Gratuiertenkolleg “Non-equilibrium phenomena and phase transitions in complex systems” and thanks the University of Crete and the Research Center of Crete for hospitality.
       
  
E.H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963). E.H. Lieb, Phys. Rev. 130, 1616 (1963). T. Tsuzuki, J. Low Temp. Phys. 4, 441 (1971). V.E. Zakharov and A.B. Sabat, Sov. Phys. JETP 37, 823 (1973). P.P. Kulish, S.V. Manakov, and L.D. Faddeev, Theor. Math. Phys. 28, 615 (1976). M. Ishikawa and H. Takayama, J. Phys. Soc. Japan 49, 1242 (1980). S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Sengstock, A. Sanpera, G.V. Shlyapnikov, and M. Lewenstein, Phys. Rev. Lett. 83, 5198 (1999). J. Denschlag, J. E. Simsarian, D. L. Feder, Charles W. Clark, L.A. Collins, J. Gubizolles, L. Deng, E.W. Hagley, K. Helmerson, W.P. Reinhardt, S.L. Rolston, B.I. Schneider, and W.D. Phillips, Science 287, 97 (2000). A.E. Muryshev, H.B. van Linden van den Heuvell, and G.V. Shlyapnikov, Phys. Rev. A 60, R2665 (1999). D.L. Feder, M.S. Pindzola, L.A. Collins, B.I. Schneider, and C.W. Clark, Phys. Rev. A 62, 053606 (2000). B.P. Anderson, P.C. Haljan, C.A. Regal, D.L. Feder, L.A.
      
Collins, C.W. Clark, and E.A. Cornell, Phys. Rev. Lett. 86, 2926 (2001). C.A. Jones and P.H. Roberts, J. Phys. A: Math. Gen. 15, 2599 (1982). N. Papanicolaou and P.N. Spathis, Nonlinearity 12, 285 (1999). F. Dalfovo and S. Stringari, Phys. Rev. A 53, 2477 (1996). E. Zaremba, Phys. Rev. A 57, 518 (1998). G.M. Kavoulakis and C.J. Pethick, Phys. Rev. A 58, 1563 (1998). S. Stringari, Phys. Rev. A 58, 2385 (1998). M.R. Andrews, D.M. Kurn, H.-J. Miesner, D.S. Durfee, C.G. Townsend, S. Inouye, and W. Ketterle, Phys. Rev. Lett. 79, 553 (1997); ibid 80, 2967 (1998). R.J. Donnelly, Quantized Vortices in Helium II, (Cambridge University Press, 1991). C.A. Jones, S.J. Putterman, and P.H. Roberts, J. Phys. A: Math. Gen. 19, 2991 (1986).