On creation and evolution of dark solitons in Bose-Einstein condensates V.A. Brazhnyi1, A.M. Kamchatnov1,2, and V.V. Konotop1 ∗ † ‡ 1Centro de F´ısica da Mat´eria Condensada, Universidade de Lisboa, Complexo Interdisciplinar, Av. Prof. Gama Pinto 2, Lisboa 1649-003, Portugal 2Institute of Spectroscopy, Russian Academy of Sciences, Troitsk 142190, Moscow Region, Russia Generation of dark solitons from large initial excitations and their evolution in a quasi-one- 3 dimensional Bose-Einstein condensate trapped by a harmonic potential is studied analytically and 0 numerically. In the case of a single deep soliton main characteristics of its motion such as a fre- 0 quency and amplitude of oscillations are calculated by means of the perturbation theory which in 2 the leading order results in a Newtonian dynamics, corrections to which are computed as well. It n is shown that long-time dynamics of a dark soliton in a generic situation deviates substantially a from outcomes of the naive application of the Ehrenfest theorem. We also consider three different J techniques of controllable creation of multi-soliton structures (soliton trains) from large initial ex- 7 citationsand calculate theirinitial parameters (depthsandvelocities) with theuseofageneralized 1 Bohr-Sommerfeld quantization rule. Multi-soliton effects are discussed. ] PACSnumbers: 03.75.Lm,03.75.Kk,05.45.Yv t f o s I. INTRODUCTION not happen in the case of bright solitons in a BEC with . t negativescatteringlength,whichinaquasi-1Dcasewere a recently observed experimentally [8]). m Experimental observations of dark solitons in a cigar- In the homogeneous NLS equation, which is exactly - shaped Bose-Einstein condensate (BEC) of sodium [1] d integrable [9], the long-time evolution of an excitation and of 87Rb [2] atoms and in a near spherical BEC of n can be reduced to the motion of a set of solitons and 23Naatoms[3]havestimulatedintensivetheoreticalstud- o to propagation of linear waves. Solitons in that case are c ies devoted to generation and evolution of such excita- well defined objects and their parameters can be found [ tions in BEC’s. At low enough temperature the conden- by means ofthe inversescatteringtransformmethod[9]. sate evolutionis described sufficiently wellby the Gross- 1 In particular, constant velocities of solitons (as well as v Pitaevskii (GP) equation [4], which is originally three- other soliton parameters) can be calculated from the as- 9 dimensional (3D) but in some cases is reducible to a 1D sociated linear spectral problem using the initial distri- 1 nonlinear Schr¨odinger (NLS) equation with an external bution of the macroscopic wave function. In the pres- 3 potential. As such cases we can mentioned a BEC with ence of a trap this method cannot be applied without 1 a “pancake” geometry (see e.g. [5]) and a low-density 0 some reservations. Moreover, a dark soliton against a BEC in a “cigar-shape” trap (see e.g. [6]). If the re- 3 nonuniformbackgroundis not a well defined object any- spective size of the condensate is much bigger than the 0 more. Therefore separation of the excitation evolution / characteristic dimension of an excitation in it, then it is into “soliton part” and “background part” is somewhat t a commonly believed that the conventional (i.e. without conditional. Even if such separation is meaningful at m potential) NLS equation can be used for understanding the initial moment of time, during propagation soliton the excitation evolution during initial intervals of time. - parameters become functions of time. In particular, a d Suchapproximationwasusedfordescriptionofdarksoli- change of a soliton shape can be so dramatic that the n tondynamicsinquasi-1DBECwithapositivescattering distinction between the soliton and the backgroundmay o length (see e.g. [7]). If, however, a size of an excitation c looseitssenseandmustbereconsidered. Thenonemeets is comparable with the size of the condensate or if one : a problem of long-time description of the soliton evolu- v is interested in long-time dynamics, then influence of a tion. This problem becomes even more complicated in i X trap potential must be taken into account. This is es- the case when an initial excitation results in formation pecially important in the case of dark solitons. One of r ofmanysolitons(solitontrains). Thus,oneoftheaimsof a the physical reasons for that is that in the presence of a the present paper is to consider evolution of excitations trap potential the BEC density becomes a function of a whichinitially canbe classifiedas one-and multi-soliton space coordinate and/or time. From the mathematical ones in 1D BEC confined by a harmonic potential. point of view it means that the trap potential changes Several aspects of this problem have already been ad- boundary conditions for the macroscopic wave function dressed in literature. In particular, there have been re- of the condensate at the infinity. (Note that this does ported several evidences that the dynamics of a dark soliton or a few solitons in a BEC trapped in a har- monic potential is oscillatory [10, 11, 12, 13]. In one of the first publications [10] it has been suggested that ∗Electronicaddress: [email protected] †Electronicaddress: [email protected] for description of the dark soliton evolution one can em- ‡Electronicaddress: [email protected] ploy the Ehrenfest theorem which allows one to define 2 a frequency of oscillations, and on this basis an empiric trary initial excitation in a nonuniform condensate con- formula mx¨ = ∂U/∂x, where U(x) is a trap poten- fined by a harmonic potential as well as their evolution. − tial, for Newtonian dynamics of a dark soliton has been The organizationof the paper is as follows. In Section proposed (although a definition for the coordinate of a IIthedynamicsofasingledarksolitoninaharmonictrap dark soliton x has not been given). A similar equation potentialisconsideredindetail. Firstofallwearguethe for a dark soliton motion against a background was de- choice of the analytical form of the background which rived in [12] by means of the multi-scale analysis. In is necessary for a stable long-time dynamics of a soliton numerical simulations provided in [13] authors observed making it nearly integrable. We show also that one can an oscillatory behavior of small amplitude dark solitons. define two characteristiccoordinates associated with the That paper however has left open a question about the dark soliton: a position of the center of mass (the mass frequency of the oscillations and their amplitude. Also, being negative)andapositionofa localminimumofthe although an important question about a choice of a par- intensity. The both values being the same in the case ticular form of the background has been posed (since a of a homogeneous background display different behavior “wrong” choice results in rapidly growing oscillations of subjecttotheeffectofthepotential. Thefirstoneisgov- the background),it has been solved by numerical means ernedbytheEhrenfesttheoremwhileanotheronecanbe on the basis of elimination of the mentioned oscillations describedon the basis ofthe perturbationtheory of soli- of the background, leaving analytical fundamentals for tons [16]. Theoreticalpredictions are comparedwith the such a choice open, as well. Moreover in Ref. [13] it has numerical simulations. Non-adiabatic effects of the soli- beenfoundthattheamplitudeofasmall-depthdarksoli- tondynamicsandthebehaviorofitsphasearediscussed. ton increases as the soliton approaches a turning point, In Section III we describe generation of soliton trains in which does not corroborate with the previous findings a trapped BEC from initial excitations of different types [10, 11, 12] where it was reported that the shape of a and study their evolution. In the case of perturbation dark soliton is preserved during propagation against a of the condensate density by a large and smooth ini- background. Finally, it is worthwhile to mention that in tial pulse we find initial parameters of created solitons Refs. [10, 11, 12, 13] one can find no information about with the use of generalized Bohr-Sommerfeld quantiza- long-time behavior of a dark soliton in a trap potential. tion rule [17]. Using the results of Section II we predict Another related issue of the theory is generation of locations of turning points which are well confirmed by dark solitons. The discussion in the literature has been directnumericalsimulations. Wealsodiscussbehaviorof mainly concerned with such methods as laser induced dark solitons generated by the phase imprinting method Ramantransitionbetweentwointernalstatesofthecon- and creation of solitons during collisions of two conden- densate [14], “phase imprinting” [2, 15], modulational satesinthepresenceoftheharmonictrappotential. The instability [6], and the use collision of two initially sepa- outcomes of the theory are summarized in Conclusion. rated condensates [10, 11]. For the sake of convenience a summary of some techni- cal results of the perturbation theory for dark solitons is InRef. [10]itwassupposedthatgenerationofthedark given in the Appendix. solitons with different initial velocities can be achieved in a collision of two pieces of initially displaced conden- sates. By numerical calculations it was shown that in II. MOTION OF A SINGLE DARK SOLITON IN the presence of a trap all solitons have essentially the A PARABOLIC TRAP same classical period of oscillations given by harmonic classical dynamics. The same mechanism of generation A. Statement of the problem ofthemulti-solitonstructureswasconsideredinRef.[11] (although a number of the generated solitons and their initial parameters were not calculated). In Ref. [14] a As it is customary, we start with the GP equation for scheme to create dark solitons in a BEC with the use of the order parameter ψ ψ(r,t): ≡ coherent Raman process which couples internal atomic ∂ψ ~2 levels of the condensate with a laser was proposed. In a i~ = ∆ψ+V (r)ψ+g ψ 2ψ, (1) trap 0 ∂t −2m | | recentworkonthe controlledgenerationof darksolitons by phase imprinting method [15] the authors provided a whereweusethestandardnotations: g =4π~2a /m,a 0 s s simple theoreticaldescription of the creationof the dark being the s-wave scattering length, which is considered solitons in the experiment [2]. Considering a very short positive,andmbeing theatomicmass;V (r)isatrap trap time scale, they neglected the trap potential, and in this potential. case for a given initial phase step, a number of solitons A self-consistent reduction of Eq. (1) to a 1D NLS andtheirinitialvelocitieswerecalculatedanalyticallyby equation can be made in various situations (see e.g. [5, mapping Zakharov-Shabat eigenvalue problem into the 6]). pendulum problem, which is mathematically much sim- (i)InthecaseofapancakeBECwesupposethatinthe pler. transverse direction (i.e. in the direction orthogonal to Thus,thesecondaimofthepresentpaperisthestudy thex-axis)thesizeofthecondensateislargeenoughtobe of generation of multi-soliton structures from an arbi- consideredinfinite in the firstapproximation. This leads 3 to the trap potential of the form V = mω2x2 where where C(ν) is a very slow function of ν equal to C(ν) trap 2 0 ≃ ω istheharmonicoscillatorfrequency. Then,inorderto 2.67 with accuracy better than 1% in the interval 0.1 0 ≤ rewrite the dynamical equation in a dimensionless form ν 0.3. Comparison of Eqs. (4) and (8) shows that ≤ we make a substitution ψ(r,t)=241(4πνasa20)−21 exp ik⊥r⊥−i~2km⊥2 t Ψ(x,t), ν ≃(cid:18)2.674·π21/4(cid:19)2/3 (asa01n0)2/3 ≃ (asa00n.40)2/3. (9) (cid:18) (cid:19) (2) where r = (y,z) and a2 = ~ , and make a change of This formula determines ν in terms of the physical pa- ⊥ 0 mω0 rameters of the system under consideration. As an ex- independent variables x ν1/2a0x/21/4, t 21/2νt/ω0, ample one can consider = 105 atoms of 87Rb (a 7→ 7→ s which results in the canonical form of the NLS equation N ≈ 5.8nm)condensed in a trap with a 1µm (what corre- 0 with a parabolic potential sponds to the frequency ω 7 102≈Hz) and with trans- 0 ∼ · verse radius 14µm. Then we get ν 0.2. As it is ∂Ψ ∂2Ψ 1 ≈ ≈ i + 2Ψ2Ψ= ν2x2Ψ (3) evident, increasing of the density by a factor 10 results ∂t ∂x2 − | | 2 in decreasing of the small parameter by approximately the same factor. (in what follows Ψ is also referred to as a macroscopic wave function). Notice that in contrast to the generally Thus, ν can be viewed as a natural small parameter accepted renormalization, here we introduced a param- of the problem and therefore in what follows we shall eter ν which on the one hand characterizes a strength consider the case ν 1. ≪ of the parabolic potential in dimensionless equation (3), In order to constructthe theory yet another condition and on the other hand is connected with the density of is to be imposed: the longitudinal dimension of the con- particles in the condensate. To elucidate this connec- densate ν−1 must be sufficiently large compared with ∼ tion, we note that the condensate wave function Ψ(x,t) a characteristic width l 1 of a typical soliton solution. ∼ is normalized according to Then it makes sense to speak about a soliton against an inhomogeneous background. More specifically, we re- ∞ |Ψ(x,t)|2dx=4π2−1/4ν1/2asa0n0, (4) qsiunicreeal/lrae0ad∼y√atν.νThi0s.1hotwheevreerlaitsionnotlasatr0o/n3g(wrehsticrhictcioorn- Z−∞ responds to the ac≈tual experimental s≈ettings) is verified. where n0 stands for the “transverse density” of particle, (ii) Considering a cigar-shaped BEC one uses the i.e. n0 = /S, isthetotalnumberofparticlesandSis multi-scale expansion (see e.g. [6] for details) with re- N N theareaofthetransversecrosssectionofthecondensate. spect to the small parameter ǫ = 8π a a /a2 1 (it (Formally one has to consider the “thermodynamical” can be viewed as smallness of the eNnergsy⊥of 0tw≪o-body limit , S at n0 =const.) interaction compared with the kinetic energy) where N →∞ →∞ Ontheotherhandletusconsideranunperturbedcon- a = (~/mν )1/2 is the linear oscillator length in the densate wave function Ψ(x,t) = exp( iµt)Ψ(x) (µ is a tr⊥ansverse dir⊥ection and ν is the linear oscillator fre- renormalized chemical potential) whi−ch solves the sta- quency in the transverse di⊥rection. The smallness of the tionary equation parameter ǫ provides the conditions when the consider- ation can be restricted by the lowest mode in the three- d2Ψ 1 +µΨ 2Ψ2Ψ= ν2x2Ψ, (5) dimensional parabolic trap: at larger ǫ interactions of dx2 − | | 2 modes become essential and must be taken into account Ψ(0)=1, Ψ′(0)=0. what leads to a set of coupled nonlinear Schr¨odinger equations. In this case the time is measured in units This is clear that Ψ decays exponentially at x , 2/ν and spatial coordinate along the cigar axis is mea- | | → ∞ and thus one can expect that sure⊥d in units a . The parameter ν is now defined as ν =a /a andc⊥anbe easily made as smallas necessary. 0 ∞ Ψ(x)2dx √2µ/ν Ψ (x)2dx= (2µ)3/2 (6) Meant⊥ime, in that case one has to verify the condition TF | | ≃ | | 3ν ǫ 1 necessaryfor applicability ofthe theory developed Z−∞ Z−√2µ/ν be≪low. Let us do that for the data analogous to ones re- where ported in [1] imposing however a lower particle density, i.e. considering a BEC of = 5 105 sodium atoms 1 N · Ψ (x)= 2µ ν2x2 (7) (a 2.7nm). Then for a 450µm and a 10µm TF 2 − onse≈estimates ν 0.1. T0hu≈s in the case at⊥h≈and the p ≈ isthecondensatewavefunctionintheThomas-Fermiap- limit of small effective parameter ν is also reachable for proximationandµ 2forν 1. Morecarefulnumerical available experimental settings. ≃ ≪ study of Eq. (5) gives corrections to (6) In order to complete the statement of the problem we have to specify the terminology “long-time” dynamics. ∞ Ψ(x)2dx= C(ν), (8) We willdo this for a pancakegeometry. A unity of time, | | ν ∆t = 1, in the dimensionless variables corresponds to Z−∞ 4 ω /ν real seconds. This means, for example, that con- achievedby requiring F(x) to be an eigenfunction of the 0 sidering a BEC of 87Rb atoms in a trap with the lon- nonlinear spectral problem [c.f. (5)] gitudinal size of order of a 1µm and assuming that 0 ∼ (νx)2 l/a0 0.3, a unity of the dimensionless time will corre- F + ω F 2ρ F3 =0, (12) ∼ xx b 0 spondto0.18ms. Atypicallife-time ofaBECreachable − 2 − (cid:18) (cid:19) intodayexperimentscanbe estimatedasatleast50ms, lim F(x)=0, (13) which means that it is meaningful (in a BEC context) x →±∞ to investigate the dynamics of the excitations for times which satisfies the following normalization conditions at least up to 300 dimensionless units. As one expects ν to be of order of the frequency of soliton oscillations in F(0)=1, Fx(0)=0. (14) a trap,the aboveestimate correspondsto approximately In (12) ω is an eigenvalue. 20 periods of oscillations. Thus, long-time dynamics in b Indeed,substitutingansatz(11)inEq.(3),oneobtains thepresentpaperisunderstoodasadynamicsdisplaying more than 10 oscillations. iΦ +Φ 2(Φ2 ρ )Φ=R[Φ,F], (15) t xx 0 − | | − B. Choice of the background and near-integrable where dynamics R[Φ,F] 2(lnF) Φ +2Φ(Φ2 ρ )(F 2 1). (16) ≡− x x | | − 0 | | − When ν =0, Eq. (3) reduces to the conventionalNLS Next,onecanestimatetheorderofmagnitudeofR[Φ,F] equation which admits a stable constant-backgroundso- considering xν small enough: lution Ψ(x,t)2 = ρ (=const), against which one can 0 constru|ct a da|rk soliton solution [9] (lnF)x ∼ν2x, Φx ∼η, (17) (Φ2 ρ )(F 2 1) ν2η2x2, 0 | | − | | − ∼ 1+eiϑexp η [x X(t)] Ψ (x,t) √ρ { 0 − }. (10) where it is taken into account that in a vicinity of the s 0 ≡ 1+exp η [x X(t)] 0 center of the potential, i.e. at xν 1, F(x) can be { − } ≪ expanded in the Taylor series Here X(t) = Vt+X0, V = 2√ρ0cos(ϑ/2) and η0 = − 2√ρ0 sin(ϑ/2) are the soliton coordinate, velocity, and F(x)=1+ 2ρ0−ωbx2+O(x4). (18) width, respectively. X0 is the coordinate at initial mo- 2 ment of time. Due to the symmetry of the problem, Substitution of(18)into (12)with subsequentexpansion the parameter ϑ which characterizes the phase differ- in powers of “small” x (i.e. x = o(ν 1)) and small ν ence of Ψ (x,t) at can be restricted to the inter- − s ±∞ yields the estimate for the eigenvalue val 0 ϑ π. For ϑ = π the BEC density vanishes ≤ ≤ at the soliton center and in this case a soliton is often ν2 called“black”while forotherϑasolitonisreferredtoas ωb =2ρ0+ +O(ν4). (19) 4ρ 0 “grey”. The limit ϑ 1 correspondsto small amplitude ≪ solitons. As it follows from the definition of the soliton width, When the trap potential is switched on, the function we have η =O(1) whenever ρ0 =O(1). This means that Ψ (x,t) given by (10) is not a solution of (3) anymore. ν 1 is indeed a small parameter of the problem in Hoswever, one could expect that if a region of soliton lo- the≪sense that R[Φ,F] = O(ν2). Then the evolution of calization is much less than the size of the background, thefunctionΦ(x,t)willbedescribedbynearlyintegrable i.e. when ν η , and if in the vicinity of the poten- equation (15) and in this sense will be sufficiently close 0 ≪ tial center the initial condition for Eq. (3) is chosen to to the evolution of Ψs(x,t) given by (10) at least for be close enough to the dark soliton one, then a solution period of time defined by t ν−4. We will refer to the ≪ of Eq. (3) can be searched in a form of a “dark soliton” respective dynamics as nonlinear. If, however, xν & 1 Φ(x,t) against an inhomogeneous background F(x), i.e. then the first term in (16) becomes exponentially small in the form (sinceasitfollowsfrom(12)decayofF(x)atx is →±∞ exponentiallyfast)andthesecondtermisapproximately Ψ(x,t)=F(x) Φ(x,t), (11) canceledwiththenonlineartermin(15). Inthatcasethe · dynamics is governed by the effectively linear equation. where Φ(x,0)=Ψ (x,0), Ψ (x,t) being given by (10). s s Ansatz (11) is meaningful if the dynamics of Φ(x,t) is close, in some sense, to the dynamics of Ψ (x,t). Since C. Nonlinear dynamics of a dark soliton. s Ψ (x,t)isgovernedbytheunperturbed(i.e. withν =0) s NLS equation, to satisfy the last requirement it is nat- Suppose that X(t) = o(ν 1). In this case estimates − ural to choose F(x) such that the resulting equation for (17) hold for the whole spatial region and one can sim- Φ(x,t) would be close to the NLS equation. This can be plify the expressions for the background with the use of 5 (18) and (19): where 1 ν2 R = (x X(t))φ x2(φ2 ρ )φ F =F0(x)+O(ν4x4), F0(x)≡1− 8ρ x2 (20) 2ρ0 − x− | | − 0 0 (cid:2) ν2f2(t) e + 4if(t)φ2φ¯ + 2φφ 2 φ¯(φ )2 and rewrite R[Φ,F] in the form x ρ | x| − x 0 ν4f3(t) (cid:0) (cid:1) ν2x i φ φ 2 (25) R[Φ,F] R[Φ,F ] Φ xΦ(Φ2 ρ ) . (21) − 2ρ2 x| x| ≈ 0 ≡ 2ρ x− | | − 0 0 (cid:21) 0 (cid:2) (cid:3) where the terms proportional to the second and third Now Eq. (15) can be treated by means of the pertur- orderoff(t)insomecasesmaybenotsmall(seebelow). bationtheoryfordarksolitons[16]. Inthisway,however, The perturbation theory can be applied to (24), (25). one meets a problem: the term proportional to xΦ in x Tothisendwepasstonewindependentvariables(x,t) the righthandside of(21)belongsto the classofpertur- → (Θ,t) where bations which effect on the soliton dynamics cannot be describedbythe adiabaticapproximationonly. Toavoid Θ=η(t)(x X(t)) and X(t)=Vt+x (t) (26) 0 this difficulty, we make an ansatz − and η(t) and x (t) are allowed to be dependent on time 0 ν2 ∂φ with the initial conditions η(0) = η and x (0) = X , Φ=φ i f(t) , (22) 0 0 0 − 2ρ ∂x and look for φ in the form 0 where φ(x,t)=eiν2ϕ(x,t) φ +ν2φ + , (27) ad 1 ··· t where (cid:0) (cid:1) f(t)= X(t′)dt′+f0, (23) Z0 φad =√ρ0eiϑ2 cosϑ +isinϑtanhΘ (28) and the constant f is to be determined later. 2 2 2 0 (cid:18) (cid:19) Necessity of the renormalization (22) can be justified istheadiabaticapproximation. Asithasbeenmentioned by full invoking the perturbation theory (see Ref. [16] above, η(t) and x (t) are (slow) functions of time, de- for details). In order to give here some simple indication 0 scribing variations of the width/depth and velocity of on the origin of the phenomenon, we notice that secu- the soliton, respectively. Quantity X(t), which gives a lar terms appearing due to perturbation of a nonlinear positionof the minimum ofthe solitonintensity must be equationareeliminatedintheso-calledadiabaticapprox- interpretedasacoordinateofthesolitoncenter. Byintro- imationwhichisobtainedfromtheexactsolitonsolution ducing the rapidly varying phase ϕ(x,t) one can satisfy by allowing parameters to be dependent on time. If in our case one substitutes e iϑ/2Ψ (x,t), where Ψ (x,t) is the equation of balance of the momentum [16] − s s givenby(10),intothelefthandsideof(3)andcomputes theimaginarypart,oneensuresthatitisanevenfunction P = 1 ∞(φ φ¯ φφ¯ )dΘ, (29) Θ Θ 2i − of the spatial coordinate. Thus if Ime−iϑ/2R[Φ,F] has Z−∞ an odd (with respect to x) component, a secular terms which reads originated by such a term cannot be eliminated by any modification of the soliton parameters. This requires in- dP = ∞(Rφ¯ +R¯φ )dΘ (30) troducingadditionalphasefactorϕ(x,t)[see(27)below]. dt − Θ Θ Thepeculiarityofsuchaphaseshift,however,isthatthe Z−∞ e e imaginary part of the respective term (after multiplying (all functions in the integrands are considered to be de- by exp( iϑ/2)) is proportional to tanhη(x X(t)) and pendent on the variables (Θ,t)). thus is −zero at x = X(t). Thus if Ime iϑ−/2R[Φ,F] is In order to find dependence of η(t) on time, one can − not zero at x = X(t), it cannot be eliminated by any use the conservation law modification of the adiabatic theory. The ansatz (22) allowsone to countthe mentioned “dangerous”term ex- dN˜ =2 ∞ dxIm(φ¯R), (31) plicitly. Taking into account the explicit form of that dt Z−∞ term, namely the fact that it is related to translational where e invariance of the system and is localized about the dark soliton kernel, it can be interpreted as an internal mode N˜ = ∞ dx(ρ φ2). (32) ofadarksolitonexcitedbyspatiallyvaryingbackground. 0 −| | Now the dynamical equation for φ with the accuracy Z−∞ O(ν4) reads Direct calculations allow one to ensure that (31), (32) give dη/dt = 0, and thus the amplitude of the soliton iφ +φ 2(φ2 ρ )φ=ν2R, (24) is conserved: η η . This results gives an explanation t xx 0 0 − | | − ≡ e 6 of the numerical findings reported in [10, 11, 12]. It is where important,however,thatithasbeenobtainedinthelimit 2V2 η ν and thus cannot be compared with the outcomes ν˜2 =ν2 1+ (38) of≫[13]. (cid:18) 9 ρ0 (cid:19) Now, computing R and substituting the result in (29) is the renormalizedfrequency and in the right hand side and(30),oneobtainssubjecttotheboundaryconditions X(t) has been substituted by X (t) given by (36). Then 0 limx ϕ(x,t)=0e: one computes | |→∞ ν2η0 Θ t 541 V2 V ϕ= X(t)dt . (33) − 4ρ0 cosh2(Θ/2)Z0 ′ ′ X(t) = (cid:18)1+ 2160 ρ0 (cid:19) ν sin(ν˜t) Let us consider a soliton trajectory. After direct sub- V3 23 1 sin(2νt)+ sin(3νt) . (39) stitution ofEq.(25)into Eqs.(A2), (A3) andthen substi- − ρ ν 216 80 0 (cid:18) (cid:19) tution of the obtained result into (A1) one can get (see Appendix A) Theobtainedresultrevealsthreeimportantfeaturesof the dynamics when the soliton velocity increases. First, dx ν2 ϑ t thefrequencyofoscillationsincreasescomparedwiththe 0 = 2sin2 +1 X(t)dt dt − 3 2 ′ ′ frequency of the harmonic trap (in the case at hand by ν2V X(cid:18)2 4ν4V sin(cid:19)2Zϑ20 tX(t′)dt′ 2 (19bVρy02t)h.eSevcaolnude,t5h4e1aVm2)p.litTudheirodf,otshceilrleateioxnisstaslsaosinmcarellas(es − 8ρ0 − 9ρ0 (cid:18)Z0 (cid:19) 2160 ρ0 ∼ ν) frequency mismatch between second harmonic and a 2ν6sin4 ϑ t 3 + 2 X(t)dt . (34) double first harmonic and between third harmonics and ′ ′ 15ρ0 (cid:18)Z0 (cid:19) the triple first harmonic. This leads to slow (compared with the period of soliton oscillations) variation of the Requiring the initial soliton velocity to be V, i.e. dx0 =0 one obtains the constant introduced in (23): amplitude of the periodic motion. Below we will observe dt |t=0 allthesephenomenainnumericalsimulations(seeFig.1). f = (6+π2) V +O(ν2V2). 0 24 ρ0η02 Below we will argue that the theory is applicable for relatively small velocities, in particular for V ν. This D. On the Ehrenfest theorem for dark solitons ∼ meansthatwearetorestricttheconsiderationtoϑclose to π. Let us first consider ϑ π ν. This allows one As was indicated in [10, 12], the frequency of the soli- | − | ∼ to verify, using Eq. (34), that X(t) in the leading order ton oscillations can be found with help of the Ehrenfest with respect to ν (designated below as X0(t)) satisfies theoremwhichmustbewrittenforthe“centerofgravity” the equation of the harmonic oscillator of the whole condensate dd2Xt20 +ν2X0(t)=0 (35) x(t)= N1 ∞ x|Ψ(x,t)|2dx, (40) Z−∞ and thus where V X (t)= sin(νt) (36) 0 ν N = ∞ Ψ(x,t)2dx (41) | | (where the initial conditions X (0)=0, X˙ (0)=V have Z−∞ 0 0 been used). Thus, the soliton center, defined as a local is a renormalized number of particles of the whole con- minimum of the intensity undergoes oscillatory motion densate (N ). Then the coordinate of the dark soli- ∼N with the frequency which is equal to the frequency of ton can be defined as the trap ν. This result also corroborates with previous investigations[10,12,13]. NoticethatfromEq. (36)and x (t)= 1 ∞ x[F(x)]2(Φ(x,t)2 ρ )dx, (42) s 0 N | | − the condition |X(t)|.1 follows the relation V ∼ν used s Z−∞ above. where Let us now return to expansion (20) and observe that IitnitshsitsilclavsaelitdheifswoleitcoonnsviedleorciXty∼isνa−llo1/w2ed(atnodbνe1/o2f≪ord1e)r. Ns = ∞(|Φ(x,t)|2−ρ0)dx (43) ofν1/2. Thendifferentiationof(34) withrespectto time Z−∞ yields (instead of (35)) can be interpreted as a negative mass of a dark soliton. One can easily show that x (t) coincides with X (t), s 0 d2X νV3 23 1 defined in subsection IIC, only in the leading order and +ν˜2X(t)= sin(2νt)+ sin(3νt) dt2 ρ 72 10 at initial stages of evolution, and thus cannot be treated 0 (cid:18) (cid:19) (37) as a coordinate of the soliton center in a general case. 7 Indeed,itisstraightforwardtoshowthatx (t)obeysthe Third, one observes that the amplitude of oscillations s exact “Newton law” are larger then V/ν predicted by Ehrenfest theorem and displaystheperiodicmodulation. Allthesedifferencesin- d2xs = ∞ ∂Vtrap Ψ(x,t)2dx, (44) creasewith growthof the initial soliton velocityV. This dt2 − ∂x | | behavior is also in qualitative agreement with the result Z−∞ (39) obtained within the framework of the perturbation which in the case at hand is reduced to the equation for theory. the harmonic oscillator. Then, taking into account the initial condition in the form Ψ(x,0)=F(x) Φs(x,0), (45) FIG.1: DarksolitontrajectoriesX(t)(solidlines): forρ0 =1 · and different initial parameters (a) ϑ = 2; (b) ϑ =2.5; (c) 0 0 which corresponds to the form of the solution given by ϑ0 = 3, and the trajectories following from the Ehrenfest (11), one computes theorem x¯s (dashed lines); Figure (d) shows soliton trajecto- ries for ϑ = 2 and different amplitudes ρ = 0.5;1;2 (solid, 0 0 V dotted and dashed line correspondingly). The parameter of xs(t)= sin(νt), (46) harmonic trap potential is ν =0.2. ν whichisanexactformula(noticethatEq.(36)isapprox- It is to be emphasized that all above effects become imate). easily visible after several periods of oscillations, while Thus, wehaveobtainedthat inthe lowestorderofν a during the firstperiod they closelyfollow the “Ehrenfest darksolitonwithasmallvelocityV ν undergoesoscil- trajectory”givenby(46)or(36). Thiscoincidenceofthe ∼ latory behavior which emerges both from the perturba- trajectories has been observed in Refs. [10, 11, 12, 13] tion theory for solitons and from the Ehrenfest theorem. where only few cycles were counted. Our numericalsim- Suchabehaviorisconfirmedbythedirectnumericalsim- ulations of first turning point also have shown that if an ulations(seeFig.1cbelow,whereν =0.2,andV 0.14). initial solitonis deep enough(what correspondsto small ≈ However,significant deviations from the Ehrenfest theo- valuesofV),thentheamplitudea =V/νfoundfrom theor rem are observed with growth of the initial dark soliton Eq.(36)andamplitudea calculatednumericallyasa num velocity [c.f. (46) and (39)]. local minimum of the intensity practically coincide with each other, which is illustrated by Table I. E. Numerical simulations of the one-soliton dynamics. TABLE I: Amplitudes of oscillations calculated theoretically atheor =V/ν and from thedirect simulation of Eq. (3) anum forρ =1anddifferentinitialvaluesoftheparameterϑ and To verify the above findings, we have done extensive 0 0 for two different trap potentials. numericalsimulationsofthedarksolitondynamicswhich is governed by the evolution equation Eq. (3) subject to ϑ ν =0.2 ν =0.1 0 the initial conditions (45). In Fig. 1 we present a dark atheor anum atheor anum solitontrajectoriesX(t),whichnumericallyweredefined 1 8.78 9.17 17.55 17.99 as the coordinate of the local minimum of the soliton, 1.25 8.11 8.33 16.22 16.59 fromthelong-timesimulations. Although,strictlyspeak- 1.5 7.32 7.49 14.63 14.96 ing, only the case depicted in Fig. 1c corresponds to the perturbation theory developed above, all the plots 1.75 6.41 6.56 12.82 13.09 display severalfeatures qualitatively coinciding with the 2 5.40 5.54 10.80 11.01 theoretical predictions. 2.25 4.31 4.42 8.62 8.77 First, the choice of the form of the background as an 2.5 3.15 3.23 6.31 6.41 eigenfunction of the problem (12) indeed allows one to 2.75 1.95 1.99 3.89 3.95 follow the dynamics up to hundreds of periods of oscilla- 3 0.71 0.73 1.42 1.44 tions, if initially soliton had large enough amplitude, i.e. one still can identify an entity called dark soliton. Second, the dynamics of such defined X(t) displays In Fig.2 we presentsnapshots of the darksoliton evo- relatively large deviations from that, which would follow lution for a region of parameters when the perturbation from the Ehrenfest theorem, when x¯ is identified with theory strictly speaking is not applicable (ϑ = 2 and s 0 the soliton coordinate. Namely, in all the cases the fre- ρ =1correspondtoV 1.08andthusV >ν). Onecan 0 ≈ quency of the soliton oscillations is higher than the fre- observe rather strong non-adiabatic effects which mani- quency,predictedfromthe Ehrenfesttheorem(see(46)). fest themselves in a significant deformation of a soliton The observed frequency change is of order of a few per- shape,namely informationofleadingandtrailingwaves cents ( 3% in Fig. 1b) which agrees with prediction before and behind the soliton in the “background” dis- ∼ (38). tribution. 8 mulas of the perturbation theory, i.e. Eqs. (10), (22), (27), and (28). Comparing the two sets of graphics one FIG. 2: The reduced shape of the initial distribution com- observesremarkablesimilarityofthemainfeatures(and, putedas|F(x)|2−|Ψ(x,t)|2(solidline)forϑ0 =2andν =0.2 what is more important, the main scales) of the dynam- at (a) t = 0; (b) t = 3; (c) t = 7.3; (d) t = 11; (e) t = 15.2; ics in the both cases. Three differences, however, are (f) t=18. Thin dashed lines show thebackground |F(x)|2. observable. First, the phase ofthe directnumericalsolu- tiondisplaysasymmetry,whichisnaturallyexplainedby the asymmetric non-adiabatic modulation of the soliton Some featuresofthe darksolitondynamicscanbe un- shape (c.f. Fig.2). Second,spatiallocations ofthe turn- derstood much better if one investigates in detail the ingpointsareslightlydifferent. Third,asymptoticvalues turning points. First of all, taking into account that ofthehydrodynamicalvelocityfarfromthesolitoncenter only a black soliton has zero velocity, it is natural to differ as well. These discrepancies are explained by the suppose that at the turning points the soliton becomes fact that the case under consideration does not satisfy black. The less is the soliton depth, the farer from the the conditions of the applicability of the adiabatic ap- center is a turning point, and hence the greater is the proximation, and thus one can expect that the last one amplitude of the oscillations. This correlation between can give only right orders of magnitude of the parame- variationsofthedepthandoftheamplitudebecomesev- ters and indicate correctly the main qualitative features ident if one compares the results depicted in Fig. 3 and (what is indeed observed in comparison of the two sets in Fig. 1. As it has been mentioned in the previous sec- of graphics). tion, these oscillations are due to mismatch between the the first and higher harmonics. In all numerical studies carried out we have observed oscillations similar to ones III. FORMATION OF SOLITONS FROM LARGE showninFig.1. Theycanbeexplainedbythe mismatch INITIAL EXCITATION 3(ν ν˜). We however did not observe pronounced os- − cillations due to difference in the second harmonics, i.e. The problem of evaluation of parameters of dark soli- due to 2(ν ν˜). We attribute this behavior to the fact − tons formed from a large initial excitation on a constant thatnon-adiabaticeffectsdue tothethirdharmonicsare background is formally solved by the inverse scatter- effectively much stronger than due to the second one. ing method [9]. In the framework of this method, the NLS equation is associated with the so-called Zakharov- Shabat linear spectral problem [9], and soliton param- eters are related with the eigenvalues of this problem calculated for a given initial condensate wave function Ψ(x,0). Iftheinitialdisturbanceislargeenough,sothat FIG.3: Timedependenceofthedepthofthedarksolitonfor the linear spectral problem possesses many eigenvalues, ρ0 =1, ν =0.2 and different initial parameters ϑ0 =2;2.5;3 then a well-knownquasi-classicalmethod can be applied (curves1,2,3 correspondingly). for their calculation. As was shown in recent paper [17], a generalized Bohr-Sommerfeld quantization rule is very convenient for this aim. To formulate this rule, it is convenient to introduce a new small parameter ε, ε 1, into Eq. (3) by means FIG.4: Behaviorofthehydrodynamicvelocityv=∂ϕtot/∂x of replacements x = x/ε, t≪= t/ε, ν = ν ε, so that the ′ ′ ′ inthevicinityofaturningpointforinitialparametersϑ =2 0 equation transforms to and ν = 0.3. ∆t is a time counted from the turning point. Left column isan outputofdirect numerical simulations and ∂Ψ ∂2Ψ 1 right one is formally obtained form the perturbation theory iε +ε2 2Ψ2Ψ= ν2x2Ψ, (47) ∂t ∂x2 − | | 2 (see thetext). where we have omitted for simplicity the primes in the Finally, an important effect accompanying changes new variables. Then the limit ε 1 corresponds to ≪ of the soliton velocity is a change of the total phase, formation of a large number of solitons from an initial ϕ =argΦ(x,t), of the solution which is especially sig- disturbance with parametersof orderof magnitudeO(1) tot nificant at the turning points. Our numerical studies of (see [17]). In framework of the inverse scattering trans- the phase evolution are shown in Fig. 4, where in the form method the NLS equation, i.e. Eq. (47) with the left column we depicted the hydrodynamical velocity v zero right hand side, is treated as a compatibility con- defined as v = ∂ϕ /∂x (see also the next section). As dition of two linear equations for auxiliary function χ, tot one can see, near the turning point hydrodynamicveloc- which we write down in the form ity changes dramatically making a flip (see Figs. 4b,c), 1 and after that soliton moves in the opposite direction ε2χxx = χ, χt = χx+ xχ, (48) A −2B B (Fig. 4d). In the right column we show the phase vari- ation during the reflection formally computed using for- (equivalent to the Zakharov-Shabat problem; see [18, 9 19]), where potential depending on the spectral parameter λ. Ac- cording the mentioned above independence of the eigen- iεΨ 2 Ψ valuesλ ontime, theycanbe calculatedwiththe useof = λ x + Ψ2 ε2 x , (49) n A − − 2Ψ | | − 2Ψ the initial distributions ρ(x,0) and v(x,0). Since we are (cid:18) (cid:19) (cid:18) (cid:19)x interested in such initial data which give rise to creation iεΨ = 2λ+ x. (50) of large number of solitons, this means that functions B Ψ ρ(x,0) and v(x,0) correspond to the problem (54) with The first equation (48) may be considered as a second a large number of eigenvalues. Then the quasi-classical orderscalarspectralproblemwitha given“potential”Ψ approach can be used for their calculation [17]. To clar- and λ playing the role of the spectral parameter. When ify this method, we have shownschematicallyin Fig.5 a Ψ(x,t) evolves according to Eq. (47) with ν = 0, the plot of the “Riemann invariant” eigenvalues λ of this problem do not change with time n t, and each eigenvalue corresponds to a soliton created λ+ = v(x,0)/2+ ρ(x,0) (55) from the initial pulse. To investigate the limit ε 1, let − us represent the condensate wave function in the≪form p which plays a role of the “quantum-mechanical poten- i x tial” for the problem (54). (The second Riemann invari- Ψ(x,t)= ρ(x,t)exp v(x,t)dx , (51) ′ ′ antλ = v(x,0)/2 ρ(x,0)canbe consideredinthe ε − (cid:18) Z (cid:19) same way.−) The Riem−ann invariant for the same distur- p p whereρ(x,t)hasameaningofthecondensatedensityand bance but with respect to uniform backgroundF(x)=1 v(x,t)isthehydrodynamicvelocity. Indeed,substitution is shown for comparison by thin solid line. We see of (51) into (47) yields the system that in both cases there is a “potential well”, but for the nonuniform backgroundcase the eigenvalues acquire 1ρ +(ρv) =0, imaginary part (“decay width”) due to tunnelling effect. 2 t x 1v +vv +ρ +ε2 (ρ2 2ρρ )/8ρ2 = 1ν2x, This means that dark solitons in confined condensate 2 t x x x− xx x −2 have a finite life-time τ determined by the imaginary (52) n (cid:0) (cid:1) partoftheeigenvalueλ ,whatcorrelateswiththeabove which for ε 0 takes the form of hydrodynamic equa- n → mentionedfactthatinthiscasethesolitonisnota“well- tions. For smooth enough functions ρ(x,t) and v(x,t), defined” object. Nevertheless, it makes sense to speak when aboutsolitonsinaconfinedcondensate,iftheirlife-times τ are much greater than the period 2π/ν of their os- ερ /ρ ρ, and εv ρ, (53) n x x ≃ | |≪ | |≪ cillations discussed in the preceding Section. It is clear what corresponds to neglecting the space derivatives in that“shallow”solitonswithsmallτndonotsurviveinthe , spectral problem (48) transforms into confined condensate, so that λn with values close to the A top of the “potential barrier” do not correspond to any ε2χ = (λ+v/2)2+ρ χ. (54) realsolitons. Onthe contrary,inthe case ofthe uniform xx − backgrounddiscussedin[17]alleigenvaluescorrespondto h i real solitons appearing eventually from the initial pulse. It is to be mentioned here that conditions (53) are noth- In multi-soliton problem there is also one more scale of ing but the conditions of applicability of the well-known time equal to time of formation of solitons from the ini- hydrodynamical approach when due to a relatively high tial pulse. For deep enough solitons it can be estimated density the two-bodyinteractionsarestrongenoughand by the order of magnitude as time necessary for solitons one can neglect the “quantum pressure” (see, e.g. [4]). with velocity V = 2λ (see below) to pass the dis- n n | | | | tance equal to the width d of the initial problem. For the problems under consideration, this time d/2λ FIG.5: SchematicplotoftheRiemanninvariantλ+ givenby ∼ | n| must be much less than the period 2π/ν. Thus, we are Eq.(55) (thick solid line); thin solid line shows the Riemann invariantfordisturbancewithrespecttouniformbackground; interestedintheeigenvaluesλn whichsatisfyinequalities dashed line shows background without disturbance. x∗ de- notes the position of localized disturbance and ∆F(x∗) = 1−F(x∗) is the change of the condensate density because 2π d τ . (56) of condensate non-uniformity. d is characteristic scale of the n ≫ ν ≫ 2λ n initialdisturbance. Thehorizontallinesofdifferentwidthin- | | dicate positions of eigenvalues λn and the width of each line characterizes the lifetime of the corresponding solitons (the Itisclearthatthereareno thesereservationsinthe case thickeris a line, thesmaller is the lifetime). of uniform background [17] where formally τn = and ∞ onecanwaitlongenoughto observeformationofsoliton Equation (54) has a formal analogy with a stationary with any value of λn. Schr¨odinger equation for a quantum-mechanical motion To calculate the real parts of the eigenvalues λ cor- n ofaparticleinthe“energy-dependent”potential,i.e. the responding to deep solitons, one can use the generalized 10 Bohr-Sommerfeld quantization rule TABLEII: Parameters ofsolitons createdfrom initial inten- sity disturbance 1 x+n 1 2 1 ε Zx−n s(cid:18)λn+ 2v(x,0)(cid:19) −ρ(x,0)dx=π(cid:18)n+ 2(cid:19), n λ(n) ϑ(n) a(tnhe)or a(nnu)m 0 0.41 2.29 2.74 2.48 n=0,1,2,...,M, (57) 1 0.65 1.72 4.35 4.23 with given initial distributions ρ(x,0) and v(x,0). We 2 0.80 1.28 5.42 5.54 suppose here that the integrand has only one maximum and x+ and x are the points where the integrandfunc- n −n tion vanishes: they depend on λ and arechosenso that n thatthe“initialcoordinates”ofsolitonscreatedfromthe the relationship (57) is satisfied (see Fig. 5). Analytical initial pulse cannot be identified exactly with X(0)= 0. form of each emerging soliton in an asymptotic region This is the reason why solitons created from one initial where it is well separated from other solitons (i.e. in the pulse do not reverse simultaneously their directions of limit t ) is expressed in terms of λ as follows: n motion even during the first period of oscillations in the →∞ confined condensate. This phenomenon is illustrated in ρ λ2 ρ(n)(x,t)=ρ 0− n , (58) Fig. 6, where the density ρ(x,t) and hydrodynamic ve- s 0− cosh2 ρ λ2 (x 2λ t)/ε locity v(x,t) of the condensate are shownas functions of 0− n − n x at the moments when two solitons in each train have hp i already reversed the direction of their propagation and aremovingtothecenter,whiletheothertwosolitonsare v(n)(x,t)=λ ρ /ρ(n) 1 . (59) s n 0 s − stillmovingfromthecenter. Thisprocessof“solitonsre- (cid:16) (cid:17) flection from the potential well” takes about 20% of the As it is clear, formulas (58), (59) and (10) with ϑ= ϑ , n whole period 2π/ν of their oscillations. η = η , and V = V , where the parameters are con- 0 n n nected by the relation λn = √ρ0cos(ϑn/2), represent the same one-soliton solution. Thus, the last formula al- FIG.6: SpacedistributionofthedensityofaBEC(a)andof lows one to find initial values of ϑ for solitons emerging itshydrodynamicvelocity(b)intheharmonictrappedpoten- n fromthedarkexcitationofthecondensatewithgivenini- tialwithν =0.3attimet=4.24withinitialexcitationtaken tialdistributionsofρ(x,0)andv(x,0)againstaconstant in the form of pulse(60) with ρ0 =1, α=0.8 and ε=0.3. uniform background. If the backgroundis not uniform but changes in space in the intervals of integrationin (57), then we can apply the same method with ρ replaced by the value of the FIG. 7: Space distribution of density of BEC (a) and its hy- 0 background density F2(x ) at the place x of the local- drodynamic velocity (b) in the harmonic trapped potential ∗ ∗ with ν = 0.3 at time t = 1 without (dotted line) and with ized initial excitation (see Fig. 5). (solid line) initial excitation that is taken as a phase step We have used this approachfor finding solitonparam- with ρ =1, α=4, κ=0.4 and ε=1. 0 eters for different types of initial excitations: (i) exci- tations of the density ρ(x) [1], (ii) excitation of the hy- (ii) For illustration of the process of solitons forma- drodynamic velocity v(x) (“phase imprinting method”) tionby the phase imprinting methodwe havechosenthe [2, 15] and (iii) collision of condensates [11]. initial conditions in the form (i)Inthefirstcaseofthedensitydisturbancetheinitial data were taken in the form α ρ(x,0)=F2(x), v(x,0)= . (61) cosh2(κx) 2 α ρ(x,0)= 1 , v(x,0)=0, (60) − cosh(x) Againthesolitonsparameterscanbecalculatedwiththe (cid:18) (cid:19) use of the Bohr-Sommerfeld quantization rule (57) and where the parameter α measures the strengthof the dis- their values correspond well to numerical simulation. In turbance. We have chosen the following values of the particular, the Bohr-Sommerfeldrule gives correct num- parameters: α = 0.8, ν = 0.2, ε = 0.3. The values ber of solitons and signs of their initial velocities. As of λn for the three deepest solitons calculated with the one can see in Fig. 7, at first the deepest soliton moves use the Bohr-Sommerfeldrule (57)areshownin TableII to the right but the hydrodynamic velocity distribution together with the corresponding values of ϑn and am- corresponds to its motion to the left. Only when the plitudes a(n) calculated for these solitonfrom (36) and local density minimum touches the x axis, the velocity theor foundfromnumericalsolutionof(3)withtheinitialdata distribution makes a flip and after that it corresponds (60). Thediscrepancyislessthan10%andiscaused,ap- to the predicteddirectionofthe solitonpropagation(see parently,by thefactthatinourcasethe createdsolitons Fig. 8). All solitons predicted by the Bohr-Sommerfeld did not reach yet the asymptotic values of their veloci- quantization rule can be observed during some time in- tiesV =2λ . Besidesthat,numericalcalculationsshow terval after their formation (see Figs. 8,9). However, in n n