Two-Dimensional Superfluid Flows in Inhomogeneous Bose-Einstein Condensates

Tw o-Dimensional Sup erfluid Flo ws in Inhomogeneous Bose-Einstein Conden sates Zheny a Y an 1 , ∗ V. V. Konotop 2 , A. V. Y ulin 2 , a nd W. M. Liu 3 1 Key L ab or atory of Mathematics Me chanization, Institute of Systems Scienc e, AMSS, Chinese A c ademy of Scienc es, Beijing 100190, Chi na 2 Centr o de F ´ ısic a T e´ oric a e Computacional an d Dep artamento de F ´ ısic a, F aculdade de Ci ˆ encias, Universidade de Lisb o a, Av enida Pr ofessor Gama Pinto 2, Lisb o a 1649-003, Portugal 3 Beijing National L ab or atory f or Condense d Matter Physics, Institute of Physics, Chi nese A c ademy of Scienc es, Beijing 100190, China W e rep ort a no vel algorithm of constructing linear and nonlinear p otentials in t h e tw o-dimensional Gross-Pitaevskii equation sub ject to given boun dary conditions, which allo w for exact analytic so- lutions. The o btained solutions represent sup erfluid flo ws in inh omogeneous Bose-Einstein conden- sates. The method is based on th e combination of t he similarit y reduction of the tw o-dimensional Gross-Pitaevskii equation to the one-d imensional nonlinear Schr¨ odinger eq uation, the latter allo wing for ex act solutions, with the conformal mapping of the given domain, where the flow is considered, to a h alf-space. The stabilit y of the obtained flow s is addressed. A num ber of stable and ph ysically relev ant examples are described. P A CS nu mbers: 05.45.Yv, 03.75.Lm, 42.65.Tg I. INTRO DUCTION Now a days o ne observes rapidly increa sing interest in studying nonlinear Schr¨ odinger (NLS) equations with in- homogeneous co efficients and, in particula r, in obtain- ing their exact analytical solutions for the physically relev a nt statements. Starting with the first results o n int egra ble inho mogeneous mo dels [1], this a ctivity re- ceived further dev elopment due to its relev ance in no n- linear optics [2 ] and in the mean-field theo r y of Bo se- Einstein condensates (BECs) [3, 4 ]. A p os sibility of con- structing exact so lutions w a s also rep orted for NLS equa- tions with inhomog eneous complex- v alued co efficients [5]. How ever, there ar e t w o main limitations of the pre sently av ailable results. First, most of them w ere obtained for one- dimensional o r quas i- one-dimensiona l sta temen ts and only a few res ults on multidimensional problems were rep orted, so far [4, 6]. Second, the most of the models allowing for construction of exact solutions were pos ed in the infinite domains. The first sugges tion of an algor ithm for constructing exa ct so lutions o f the one- dimens ional (1D) NLS equation on a half-line , mo deling a BEC inter- acting with a rigid surface, was rece ntly rep orted in [7]. Here we sho w that the requirements of one- dimensionality and un boundness of the do ma in ca n be remov ed and exact a na lytical solutions can be obta ined for mo dels defined o n b ounded 2D domains and de- scrib ed b y the NLS equation with inhomog eneous linear, V ext ( r ), and nonlinear, g ( r ), p otentials (th e b o th being real-v alued functions of the s patial co o rdinates). More- ov er, s o me of the rep o rted solutions ar e found to b e sta- ble, and thus having particular physical relev a nce. The paper is org anized as follows. In Sec . II, w e present the 2D physical mo del a nd giv e the conformal ∗ zyy an math@y ahoo.com mapping to reduce the 2 D ph ysical mo del s ub ject to given b oundar y conditions to nonlinear ordinary diff er- ent ial equation so lved. In Sec. I I I, we concentrate on the tw o simples t repr esentativ e exa mples to illustrate our nov el metho d. Sec. IV is dev oted to the numerical sim- ulations for the solutions obtained in the do main D 2 . In Sec. V, w e study the g eneralizatio n of the c o nformal mapping and present exact solutions o f the 2D physi- cal mo del. Finally , the outcomes are summarize d in the Conclusion. II. THE PHYSICAL MODEL AND CONFOR MAL MAPPING T o b e specific, we deal with the 2D no nlinear physical mo del i∂ t Ψ( r , t ) =  − 1 2 ∇ 2 + V ext ( r ) + g ( r ) | Ψ( r , t ) | 2  Ψ( r , t ) , (1) where r ≡ ( x, y ) ∈ D ⊂ R 2 , D is an o pe n domain, and ∇ ≡ ( ∂ x , ∂ y ). W e are particularly in terested in applica- tions of o ur results to BE C flows, where Ψ( r , t ) is the macrosco pic w av efunction and model (1) is also ter med, the Gro ss-Pitaev skii (GP ) equation [8]. W e explo re the flexibility o f p otentials in the BE C applications, i.e., p os- sibility of manipula ting them by external electric and/or magnetic fields for the sake of creation of desira ble spa tial configuratio ns fo r the linear p otential and for the scatter- ing length of th e tw o -b o dy interactions (the latter per- formed thr ough the F eshbac h reso nance technique [9]). W e also notice that the mo del (1) has also direct r ele- v ance to the mean-field theory o f exciton-p olariton con- densates [10]. There on the one hand, the 2 D sta temen t, i.e. the statement considered in this pap er, is the mos t t ypical one. On the other hand applying the exter nal pump from the free edges of a specimen one can create 2 different kinds of nonzero conditions (nonzero curr ents, as requir ed b ellow in the pr esent pa p er). W e concen trate on a B EC in a do main D bo unded b y impenetr able walls. Respectively E q. (1) will b e supplied by the zero conditions given at the b oundar y of the do- main D , which we denote a s ∂ D , i.e., w e imp ose Ψ( r ) = 0 for all r ∈ ∂ D . Our go al is to find an algo rithm allowing for sy s tem- atic co nstructions of the linear , V ext ( r ), and nonlinear , g ( r ), potentials, for which the formulated Diric hlet prob- lem allows for exa ct analytical solutions. T o this aim we assume that there exists a complex analytic function ζ ( r ) ≡ η ( r ) + i ϕ ( r ) = f ( z ) (2) of the complex v ariable z = x + iy ∈ C , whic h provides the conformal mapping of the contour ∂ D to the imaginary axis of ζ ( r ), i.e., to η ( r ) = 0, such that the domain D is mapp ed in to the right ha lf-plane in terms of the new v aria bles ( η , ϕ ): η ( r ) > 0. D ue to a nalyticity of the mapping ζ ( r ) the Ca uch y- Riemann eq uations on the pair of re al-v a lued functions η ( r ) and ϕ ( r ) hold: ∂ x η ( r ) = ∂ y ϕ ( r ) , ∂ y η ( r ) = − ∂ x ϕ ( r ) . (3) Resp ectively , the following constr aints ∇ 2 η ( r ) = 0 , ∇ 2 ϕ ( r ) = 0 , ∇ η ( r ) · ∇ ϕ ( r ) = 0 (4) are v erified, as w ell. In other w o rds, η ( r ) and ϕ ( r ) are bo th 2D ha rmonic functions with orthogo nal gr adients. Moreov er, o n the basis of the Cauchy-Riemann equa tions we find the relation |∇ η ( r ) | 2 ≡ |∇ ϕ ( r ) | 2 . W e r estrict the co ns ideration to the linea r and nonlin- ear p otentials which can be represented in terms of η ( r ) as follows V ext ( r ) ≡ − ε 2 |∇ η ( r ) | 2 , g ( r ) ≡ G 2 |∇ η ( r ) | 2 (5) where ε and G are rea l para meters. Without loss o f generality , we choose G = ± 1 . The parameter ε de- termines the pro p o rtionality co efficient betw een the p o- ten tials: V ext ( r ) /g ( r ) = − ε/ G . Then the change of the depe ndent v aria ble s Ψ( r , t ) → ψ ( η , ϕ, t ), allows one to reduce Eq. (1) to the 2D form i∂ t ψ =  ( η, ϕ )  − ∂ 2 η − ∂ 2 ϕ − ε + G | ψ | 2  ψ . (6) Here  ( η , ϕ ) is the positive definite function, defined as  ( η , ϕ ) = |∇ η | 2 / 2 ≡ |∇ ϕ | 2 / 2, where the gradients must be expressed in terms o f η and ϕ (see the examples be- low). The obtained equation (6) is co nsidered for η > 0 and has to b e supplied with the zero b oundar y condition ψ ( η = 0 , ϕ, t ) = 0. First, we concentrate on time-indep endent solutions of E q . (1): Ψ( r , t ) ≡ ψ ( η , ϕ ), wher e ψ ( η , ϕ ) s olves the 2D stationary GP eq uation with constan t co efficients εψ = − ( ∂ 2 η + ∂ 2 ϕ ) ψ + G | ψ | 2 ψ . Particular solutions of this equation can be r epresented as ψ ( η , ϕ ) = e iν ϕ φ ( η ) , (7) where ν is a constant and the real-v alued function φ ( η ) solves the problem E φ ( η ) = − ∂ 2 η φ ( η ) + G | φ ( η ) | 2 φ ( η ) , φ (0) = 0 , (8) with E = ε − ν 2 and η > 0. T urning to the ph ysical meaning of the obtained so - lutions, we observe tha t it follows from the ansatz, ψ ( η , ϕ ) = e iν ϕ φ ( η ), that ∇ ϕ ( r ) can b e iden tified as the sup e rfluid velocity . Hence, the introduced a nalytic func- tion f ( z ) is no thing but the co mplex potential of the resp ective tw o-dimensional flow. As it is well known [12] such a p otential defines the cur rent, J C of the fluid through a given contour C ⊂ D , as well as the circu- lation Γ C along C : Z C f ′ ( z ) dz = Γ C + iJ C (9) (the prime stands for the der iv ative with resp ect to z ). Since in our case f ( z ) is analytic, this in tegral is zero for any closed con tour C b o unding a simply connected domain. Th us the describ ed flow has ne ither sourc e s nor vorticit y in D . W e als o obs erve that if the change of v aria bles im- plies growth of |∇ ϕ ( r ) | with r , then the ph ysical mean- ing might hav e only s olutions with densities decaying at the infinit y (thus ens uring decaying curr ents). Whenev er one concerns with finite densities at the infinity , the phys- ically meaningful solutions w ould corr e s p o nd to ν = 0. This last constraint is assumed in what follows. Finally , we notice that if the contour ∂ D included in the pro p osed scheme is c lo sed, this implies that linea r and (or) no nlinear p otentials are div ergent at some point(s) of the bo undary . This is c le ar from the nature of the conformal mapping, since in this case there should b e a po int of the b oundar y , which is mapped into the infinity po int. In its turn such a p oint gives a n orig in to the singularity of the line a r and nonlinear po tent ials. Suc h cases will b e exclude in what follo ws, although they still may ha ve physical relev ance. II I. EXAMPLES OF EXA CT SOLUTIONS While lar ge diversit y of the domains can b e c o nsidered, here we co ncentrate on the t w o simplest repres entativ e examples for the conformal mapping (2). A. The quadran t x > | y | with the b oundary x = ± y The firs t domain considered below, is given b y D 1 = { x > | y |} (i.e., D 1 is the quadr a nt of the ( x, y ) − plane) with the b oundary ∂ D 1 = { x = y , y > 0 } ∪ { x = − y , y < 0 } . Resp ectively , the conformal mapping is chosen as f ( z ) = z 2 , and th us η ( r ) = x 2 − y 2 and ϕ ( r ) = 2 xy . The linear and nonlinear po tent ials a re now given by V ext ( r ) = − 2 ε | r | 2 , g ( r ) = 2 G | r | 2 , (10) 3 i.e., they ar e the linear expulsive parab olic p otential and the par a b olic nonlinearity . No w  ( η , ϕ ) = 2 p η 2 + ϕ 2 . B. The strip 0 < y < π with the b oundary y = 0 , π Another domain explor ed below is a strip D 2 = { x ∈ R , 0 < y < π } with the bo undary ∂ D 2 = { x ∈ R , y = 0 } ∪ { x ∈ R , y = π } . Now the function perfor ming co n- formal mapping to the upp e r half-plane is f ( z ) = e z and the new v ariables are deter mined as η ( r ) = e x sin y and ϕ ( r ) = e x cos y . The linear and nonlinear potentials al- lowing for the exa ct solutions a re now giv en by V ext ( r ) = − 2 εe 2 x , g ( r ) = 2 G e 2 x , (11) and  ( η , ϕ ) = 1 2  η 2 + ϕ 2  . C. Exact solutions T urning to the e xact solutions, below we consider only the case of repulsive interactions G = 1, as the mo s t natural candidate to pro duce stable stationary flo ws. As the t wo simplest so lutions o f Eq. (1) (see e.g., Refs. [6 , 11]) we study , the “dark soliton” shape Ψ ds ( r ) = ψ ds ( η , ϕ ) = √ E tanh r E 2 η ( r ) ! e iν ϕ ( r ) (12) and the no nlinear p erio dic mo dulation Ψ sn ( r ) = ψ sn ( η , ϕ ) = k √ 2 E √ 1 + k 2 sn √ E η ( r ) √ 1 + k 2 , k ! e iν ϕ ( r ) (13) with E = ε − ν 2 > 0 and k ∈ ( 0 , 1] b eing the mo dulus of the Jacobi elliptic s n–function [notice that ψ ds (0 , ϕ ) = ψ sn (0 , ϕ ) = 0, in conformity with the imp osed b oundary conditions]. These exact solutions g iven by Eqs. (12) and (13) and the co rresp onding velocity fields in the doma ins D 1 and D 2 are illustrated in Fig. 1. The obtained exact solutions, how ever leav e several op en ques tions r elated to their pr actical feasibilit y . First, the stability of the flo ws was not in vestigated, so far. Second, in all considered cases infinitely growing p oten- tials were used, while any cut-off (which exis ts in the real world) may strongly perturb, and even destroy the solutions. T o address these issues we now turn to direct nu merical sim ula tions. In the case at hand the p o tentials V ext ( r ) and g ( r ) grow with x , while the density go es asymptotica lly to (or is b ounded by) a certain constant lev el. Bea ring this in mind we constr uct a ph ysical system with b ounded po tent ial in the following w ay . F or x < 0 we consider shifted linear and nonlinear potentials V ext ( r − r 0 ) and g ( r − r 0 ), where r 0 = ( x 0 , 0) w ith x 0 being a constant shift vector, while V ext ( r ) a nd g ( r ) are given b y the analytical formulas (5). T o define the ph ysical potentials at x > 0 we mirror them at x = 0, thus o btaining V ext ( − r + r 0 ) and FIG. 1. (color online). The p rofi les for exact solutions given by Eqs. (12 ) and (13 ) and the correspon d ing v elocity fi elds in the domains D 1 and D 2 . Left column: Panels (a), (b), (c) show the ve locity field ∇ ϕ ( r ) = (2 y , 2 x ), and densities | Ψ ds ( r ) | 2 and | Ψ sn ( r ) | 2 for solutions given by (12) and (13) in the domain D 1 . Right column : Panels (d), (e), (f ) sho w the velocity field ∇ ϕ ( r ) = ( e x cos y , e x sin y ) , and den sities | Ψ ds ( r ) | 2 and | Ψ sn ( r ) | 2 for solutions given by (12) and (13) in the domain D 2 . The parameters are ε = 2, ν = 0, a nd k = 0 . 8. g ( − r + r 0 ). T o confine the condensate along y -c o ordinate we in tro duce additional linear tarp p otential as follows: v add = 0 for r ∈ D and v add = V 0 , whe r e V 0 is large enough, for r / ∈ D . F or more de ta ils s ee the next c hapter where we discuss the numerical studies of the problem. IV. DIRECT NUMERICAL SIMULA TIONS Now w e can p erform direct numerical s imu lations o f Eq. (1) with the additional confining po tential v add tak- ing the initial distribution of the field in the form given by the ana lytical for m ulas ψ ( r − r 0 ) for x < 0 and ψ ( r 0 − r ) for x > 0. As it is clear our ansatz do e s not sa tisfy Eq. (1) only along the line x = 0 and in the area s y < 0 and y > π . Since, how ever we imp ose strong confining po tent ial, we exp ect that the field is very w e ak outside the strip e 0 < y < π and our ansatz stays sufficiently close to the real stationary solution of Eq. (1 ) with the int ro duced physical p otentials. The a nalytically found 4 x g V ext -400 0 g 0 200 -10 0 10 0 4 x y 0 8 V ext 1 2 x 2 10 x=0 y=0 FIG. 2. (color online). The ansatz produ ced from the so lu- tions the with parameters ε = 2, ν = 0. I n the upp er part of the figure the distribution of the linear potential along y is shown for x = 0. On the righ t the distributions of the linear and nonlinear p otentials are shown alongside with th e introduced loss fo r y = 0. In the areas mark ed by “1” and “2” the anzatz exactly coincides with the analytical so lution (neglecting exp onentiall y w eak loss). solution and the corr e s po nding p o tent ials are shown in Fig. 2. Within the ar eas “1” and “2” the ansatzes used in the numerics a s the initial condition coincides w ith the analytical solutions. How ever on the bounda r y betw een these a reas the ansatzes do not sa tisfy the equation for stationary fields. The a nsatz do es not satisfy the equa - tion outside the areas either. How ever the densit y of the condensate is lo w outside the areas and so ther e are so me reasons to be lie ve that the ans atz is close to the station- ary solution. Let us remind here that in o ur n umer ical simulations to k eep the condensate lo calized within the strip e we used strong linear p otential, see the distribution of the linear potential alo ng y . W e carr ied out n umerical sim ulations fo r the so lutio ns obtained in the doma in D 2 and found out that some non-stationar y excitations appe a r but the solution sur- vive and s tay rather c lo se to the initial fie ld distribu- tion. T o c heck the stability of the s olution for longer time w e had to g et rid of the pr opagating excitations. Since w e p erfor m numerical simulations in finite win- dows the only w ay to eliminate th e propagating excita- tions is to int ro duce losses in the area where initially the densit y of the condensa te is exp o ne ntially weak. T o do this w e in tro duce linear losse s in the form γ = γ 0  exp( − ( x − x l ) 2 /w 2 0 ) + exp( − ( x + x l ) 2 ) /w 2 0  . Then the losses are negligible in the area where our analyti- cal solution is big and so the solution is practically un- affected by the artificial losses used in the numerics. In the same time all the perturba tions propaga ting aw ay ar e quickly abso rb ed. In this case the excitations disappea r FIG. 3 . ( color online). Pa nels (a) and (d) show the initial distribution of th e den sit y of the condensates for cases 1 and 2; panels (b) and (e) sho w the distributions of the condensate densities at t = 100; panels (c)and (f ) show the differences b etw een t h e initial density distributions and the density dis- tributions at t = 100. The shifts of the initial conditions and the potentials are x 0 = − 2 . 5 for (a)-(c) and x 0 = − 2 for (d)- (f ). The sol ution parameters are ε = 2, ν = 0. The field is k ept within the stripe by strong rep elling linear p otential v a = 10 4 for y < 0 and y > π . The parameters for the linear losses are γ 0 = 5, w 0 = 1, x l = 10. very quickly a nd one can see in Fig. 3 that the surviv- ing solution is very close to the initial distribution. As it is e x pe c ted the sta tio nary s o lution deviates from the ansatz most stro ngly ar ound the p oints ( x = 0, y = 0) and ( x = 0, y = π ). In the case o f attr a ctive interactions (o f a negative scat- tering length) the solutio n app ears to b e unstable aga inst collapse. V. THE GENERALIZED CONF ORMAL MAPPING Next we a ddress po ssibilities of getting mo re general t yp es of the potentials, allowing for exact solutions. T o this end we introduce the generalized relations [c.f. the Cauch y-Riemann equations (3)]: ρ 2 ( ˆ η ) ∂ x ˆ η ( r ) = ∂ y ˆ ϕ ( r ) , ρ 2 ( ˆ η ) ∂ y ˆ η ( r ) = − ∂ x ˆ ϕ ( r ) , (14) where ρ 2 ( ˆ η ) is a po sitive-definite function of ˆ η ( r ) only [as it is clear the arguments presented below can be also applied to the case where ρ 2 ( ˆ ϕ ) is a function of ˆ ϕ ( r ) only]. These equations still define transformation to the orthogo nal co o r dinates ( ˆ η , ˆ ϕ ), ho wev er now satisfying the 5 relations [c.f. Eq. (4 )] ∇ · [ ρ 2 ( ˆ η ) ∇ ˆ η ] = 0 , ∇ · [ ρ 2 ( ˆ η ) ∇ ˆ ϕ ] = 0 , ∇ ˆ η · ∇ ˆ ϕ = 0 . (15) Moreov er, now we have that |∇ ˆ ϕ ( r ) | 2 ≡ ρ 4 ( ˆ η ) |∇ ˆ η ( r ) | 2 . Notice that the generalized case [c.f. Eq. (15)] c an b e reduced to the case mentioned above [c.f. Eq. (4)] in the sp ecial case ρ 2 ( ˆ η ) ≡ 1. Next we rep eat the steps describ ed abov e for the con- formal mapping, and by the direct a lgebra show that the generalized ansa tz ψ ( ˆ η , ˆ ϕ ) = ρ ( ˆ η ) e iν ˆ ϕ ( r ) φ [ ˆ η ( r )] (16) solves the 2D stationar y GP eq uation (1) provided that φ ( r ) satisfy the s ta tionary equation (8) with η ( r ) substi- tuted b y ˆ η ( r ) and the linear and nonlinear p o tentials ar e given b y V ext ( r ) ≡ ∇ 2 ρ 2 ρ + ν 2 (1 − ρ 4 ) − ε 2 |∇ ˆ η | 2 , g ( r ) ≡ G |∇ ˆ η | 2 2 ρ 2 . (17) Notice that no w the linear and nonlinear p otentials are not pro p ortional to ea ch o ther any more [c.f. Eq. (5)]. T o construct a par ticular example, we choose ˆ η ( r ) = 2 η ( r )+ η 2 ( r ) 2 , ˆ ϕ ( r ) = ϕ ( r ) , ρ ( ˆ η ) = 1 p 2 + η ( r ) (18) with η ( r ) and ϕ ( r ) solving Eq. (4). One can ensure that this choice also satisfies the co nditions (15). Then the simplest solutio n of Eq. (8) with the r e puls ive nonlin- earity ( G = 1) is giv en by we s tudy , the “dark soliton” shap e (1 2) with η ( r ) s ubs tituted by ˆ η ( r ) and v alid for the p ositive chemical potential E = ε − ν 2 > 0. Then, according to the genera lized ansatz with ˆ η ( r ), ˆ ϕ ( r ), and ρ ( ˆ η ) given by Eq. (18), as w ell as η ( r ) a nd ϕ ( r ) defined in the domain D 2 , we obtain the “ dark soliton” ˆ ψ ds ( r ) = √ E ρ ( ˆ η ) tanh r E 2 ˆ η ( r ) ! exp( iν e x cos y ) , (19) with η ≡ η ( r ) = e x sin y , which solves Eq. (1) with the linear and nonlinear po tent ials g iven by [c.f. Eq. (17)] V ext ( r ) = e 2 x 8(2 + η ) 2  3 − 16(4 E + ν 2 ) − 2 η  ν 2 (3 + 2 η ) +2 E (4 + η )( η 2 + 4 η + 8)   , g ( r ) = 1 2 e 2 x (2 + η ) 3 . The presented s o lutions found in the domain D 2 be- hav e in a way similar to the previo us example, see Fig. 3. The numerical sim ulations prov e that the found analyti- cal solutions can indeed be used as goo d approximation for the stable stationar y s o lutions of E q. (1) descr ibing systems with ph y sically re lev ant p o tentials. VI. CONCLUSION T o conclude, we hav e shown that exa c t analy tical s olu- tions can be obtained in a lar ge clas s of t w o-dimensiona l Gross-P itaevskii equations with inhomogeneous linear and nonlinear p otentials, defined on bounded domains. The metho d o f constructing the mo dels is based on the prop erly defined co nformal mapping of the given domain int o a co mplex half-plane. In the cont ext o f applica- tions to Bo se-Einstein condensates, the obtained s o lu- tions having non trivial phase dep ending on spatial co or- dinates can be in terpreted as s upe r fluid flows. In the case of nega tive scattering leng th (repulsiv e interactions) the background flows, i.e., ones ha ving no zer os in the op en spatial doma in, a pp e ar to be stable. The obta ined re- sults genera lize previous studies devoted to construc tio n of the ex act solutions, using the s e lf- s imilar transforma- tion, to the t wo-dimensional mo dels given on bounded domains. Moreov er, the ideas pr esented in this paper are also able to apply in t wo-dimensional cubic-quint ic mo dels, tw o-dimensional m ulti-compo nent mo dels, etc., and to design linear and no nlinea r potentials for con trol of Bose-Einstein condensates and nonlinear optical fibe rs in limited spatial domains. ACKNO WLEDGMENTS Y an was supp or ted b y the NSF C under Grant No. 11071 242. VVK and A VY were supp orted by the grant PEst-OE /FIS/UI061 8/20 11 a nd b y 7th Europ ea n Com- m unity F ramework Pr ogra mme under the gr ant PIIF- GA-2009- 23609 9 (NOMA TOS). Liu was suppo rted by the NKBRSFC under Gra nt No. 2011CB9 2150 2. [1] H.-H. Chen and C.-S . Liu, Phys. Rev. Lett. 37 , 693 (1976); M . Brusc hi, D. Levi, and O. Ragnisco, Il Nuov o Cimen to A 53 , 21 (19 79); R. Scharf and A. R . Bishop, Phys. Rev. A 43 , 6535 (1991); V. V. Konotop, O. A. Chub yk alo, and L. V ´ azquez, Phys. Rev. E 48 , 563 (1993); L. Gagnon and P . Winternitz, J. Phys. A 26 , 7061 (1993); V. V. Konotop, Theor. Math. Ph ys. 99 , 687 (1994). [2] V. N. Serkin, and A . Hasega wa, Ph ys. Rev. Lett. 85 , 4502 (2000); B. A. Malomed et al ., J. Opt. B: Q uantum Semiclass. Opt. 7 , R53 (2005); S. A. Ponomare nko and G. P . Agraw al, Ph ys. Rev. Lett. 97 , 013901 (2006); V. N. Serkin, A . Hasegaw a, and T. L. Bely aev a, Phys. Rev. 6 Lett. 98 , 074102 (2007); Z. Y . Y an, Ph ys. Lett. A 374 , 672 (2010). [3] Z. X . Li ang, Z. D. Zhang, and W. M. Liu, Ph ys. R ev. Lett. 94 , 05040 2 (2005); J. Belmon te-Beitia et al., Ph ys. Rev. Lett. 98 , 064102 (2007); J. Belmonte-Beitia et al., Phys. Rev. Lett. 100 , 164102 (2008); A. T. Av elar, D . Bazeia, and W. B. Ca rdoso, Ph ys. Rev. E 79 , 025602(R) (2009). [4] V. M. P´ erez-Garc ´ ıa, P . J. T orres, and V. V. Konotop, Physica D 221 , 31 (2006). [5] F. Kh . Abdu llaev et al ., Phys. Rev. E 82 , 056606 (2010). [6] Z. Y . Y an and V. V. Konotop, Phys. Rev. E 80 , 036607 (2009); Z. Y. Y an and C. Hang, Ph ys. Rev. A 80 , 063626 (2009); S. H . Chen and J. M. Dudley , Phys. R ev. Lett. 102 , 233903 (2009); Z. Y. Y an, V. V. Konotop, and N. Akhmediev, Ph ys. Rev . E 82 , 036610 (2010); Z. Y. Y an, Phys. Scr. 78 , 035 001 (2008); Z. Y. Y an, Phys. Lett . A 374 , 4838 (2010); Z. Y . Y an, K. W. Cho w, and B. A . Malomed, Chaos, Solitons & F ractals, 42 , 3013 (2009). [7] Y u. V . Bludov, Z. Y. Y an, and V. V. Konotop, Ph ys. Rev. A 81 , 063610 (2010). [8] F. D alfov o et al ., R ev. Mo d. Phys. 71 , 463 (1999). [9] P . O. F edichev et al ., Phys. R ev. Lett. 77 , 2913 ( 1996). [10] A. Amo, D. Sanvitto, F. P . Laussy , D. Ballarini1, E. del V alle, M. D. Martin1, A. Lema ˆ ıtre, J. Bloch, D. N. Krizhanovskii, M. S. S kol nick, C. T ejedor, and L. Vi˜ na Nature, 457 , 291 (2009); A. Amo, T. C. H. Liew, C. Adrados, R. Houdr´ e , E. Giacobino, A. V. K a vokin and A. Bramati N ature Ph otonics, 4 , 361 (2010). [11] L. D. F addeev and L. A. T ak hta jan, Hamiltonian Meth- o ds in the T he ory of Solitons (Springer-V erlag, Berlin, 1987). [12] se e e.g. A. J. Ma jda and A. L. Bertozzi, V orticity and Inc ompr essible Flow (Cambridge Un ivers ity Press, 2002).