Nonlinear Schrödinger equation with spatiotemporal perturbations

Nonlinear Schrödinger equation with spatiotemporal perturbations Franz G. Mertens Physikalisches Institut, Universität Bayreuth, 95440 Bayreuth, Germany Niurka R. Quintero Departamento de Física Aplicada 1, E.U.P. Universidad de Sevilla, Virgen de África 7, 41011 Sevilla, Spain A. R. Bishop Theoretical Division and Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA 共Received 28 July 2009; published 20 January 2010兲 We investigate the dynamics of solitons of the cubic nonlinear Schrödinger equation 共NLSE兲with the following perturbations: nonparametric spatiotemporal driving of the form f共x,t兲=aexp关iK共t兲x兴, damping, and a linear term which serves to stabilize the driven soliton. Using the time evolution of norm, momentum and energy, or, alternatively, a Lagrangian approach, we develop a collective-coordinate-theory which yields a set of ordinary differential equations 共ODEs兲for our four collective coordinates. These ODEs are solved analytically and numerically for the case of a constant, spatially periodic force f共x兲. The soliton position exhibits oscillations around a mean trajectory with constant velocity. This means that the soliton performs, on the average, a unidirectional motion although the spatial average of the force vanishes. The amplitude of the oscillations is much smaller than the period of f共x兲. In order to find out for which regions the above solutions are stable, we calculate the time evolution of the soliton momentum P共t兲and the soliton velocity V共t兲: This is a parameter representation of a curve P共V兲which is visited by the soliton while time evolves. Our conjecture is that the soliton becomes unstable, if this curve has a branch with negative slope. This conjecture is fully confirmed by our simulations for the perturbed NLSE. Moreover, this curve also yields a good estimate for the soliton lifetime: the soliton lives longer, the shorter the branch with negative slope is. DOI: 10.1103/PhysRevE.81.016608 PACS number共s兲: 05.45.Yv I. INTRODUCTION The nonlinear Schrödinger equation 共NLSE兲is one of the paradigms of soliton physics, because it represents a completely integrable system and has very many applications in practically all fields of physics, which are listed and discussed in several review articles关1–3兴. For applications it is important to study the perturbed NLSE iut+uxx +2兩u兩2u=R关u共x,t兲;x,t兴.共1兲 Many different kinds of perturbations Rhave been considered and in particular the dynamics of a single soliton under these perturbations was investigated 关1,2兴. In nonlinear optics u共x,t兲is proportional to the amplitude of the electric field, i.e., 兩u共x,t兲兩⬃兩E共Z,T兲兩. Here the evolution in time is replaced by the propagation in Zand x=Tis the so-called transverse distance 关2,4兴. In real transmission systems Rmay describe the fiber loss, the gain due to the periodic amplification, the presence of filters in the fibers, etc. In this paper we consider the following combination of perturbations: R=f共x,t兲−i ␤ u共x,t兲− ␦ u共x,t兲共2兲 with the real parameters ␤ and ␦ and the nonparametric spatiotemporal driving force f共x,t兲=aeiK共t兲x,共3兲 which yields several interesting effects, as we will see. The literature so far has mostly dealt with parametric driving 关1,2,5–9兴. Nonparametric 共external兲driving was studied without space dependence, e. g., f= ⑀ exp共i ␻ t兲关10–12兴, and with a periodic space dependence, e.g., f= ⑀ exp关i共kx− ␻ t兲兴 关13,14兴. Moreover, f= ⑀ exp关ig共x,t兲−i ␻ t兴, where gis a function of x−vt, was considered, but no localized solutions were discussed 关14兴. Our driving term in Eq. 共3兲was already used in the discrete form fn共t兲=aexp关in ␾ 共t兲兴. Here ndenotes the nth resonator in an array of coupled nonlinear optical waveguides, in which discrete cavity solitons can be excited 关15兴. The array can be modeled by a discrete nonlinear Schrödinger equation 共DNLSE兲, where ␾ is the incident angle of a laser pump light 关15兴. In order to obtain a ratchet effect, a biharmonic ␾ 共t兲 was considered which breaks a temporal symmetry 关16兴. In this paper we work with an arbitrary function K共t兲in the driving term Eq. 共3兲and develop a collective-coordinate 共CC兲theory for the soliton dynamics which results in a set of nonlinear coupled ordinary differential equations 共ODEs兲for the CCs 共Secs. II and III兲. In order to obtain analytical and numerical solutions we then consider the case of temporally constant, spatially periodic driving f共x兲=aexp共iKx兲with constant K共Sec. IV兲. Although the spatial average of f共x兲 vanishes, there is transport: the soliton performs a unidirectional motion on the average, in contrast to the case of the perturbation R=V共x兲u共x,t兲with a periodic potential V共x兲in which the soliton performs an oscillatory motion around a minimum of V共x兲关7兴. Solutions of the CC equations for the case of a harmonic or biharmonic time dependence of K共t兲 will be presented in a second paper. The second term, −i ␤ u共x,t兲with ␤ ⬎0, in the perturbation 共2兲is a damping term which allows us to obtain a bal- PHYSICAL REVIEW E 81, 016608 共2010兲 1539-3755/2010/81共1兲/016608共11兲©2010 The American Physical Society016608-1 ance between the energy input from the driving and the dissipation. Other more complicated damping terms have been considered in 关17,18兴. The parameter ␦ in the third term on the rhs of Eq. 共2兲is the cavity detuning parameter in the above mentioned DNLSE 关15,16兴. In our paper, ␦ will turn out to be decisive for the stability of the driven soliton. For ␦ ⱖ0 the soliton radiates phonons 共i.e., linear excitations兲and eventually vanishes, or even breaks up into several solitons. For ␦ ⬍0 the situation is more complicated and will be discussed below. Section Vpresents a stable and an unstable stationary solution for the case without damping. For the region around the stable solution the CC theory yields solutions in which all CCs exhibit oscillations with the same intrinsic frequency ⍀. However, tests by simulations, i.e., numerical solutions of the perturbed NLSE, reveal that the oscillatory solutions are stable only for certain regions of the initial conditions. These regions become broader when ␦ is more negative 共Sec. VI兲. We conjecture that the stability of any of these oscillatory solutions can be predicted by our CC theory by calculating the curve P共V兲, where P共t兲and V共t兲are the momentum and the velocity of the soliton, respectively. This means that every point on this curve is visited during one period of the oscillatory solution. The conjecture is that the soliton will become unstable in a simulation, if the curve P共V兲has a branch with negative slope, and this is confirmed by our simulations 共Sec. VI兲. Interestingly, the curve P共V兲not only predicts whether the soliton is unstable, but it also allows us to estimate the soliton lifetime: this time is longer, the shorter the branch with negative slope is 共Sec. VI兲. The stability criteria for the NLS solitons available in the literature cannot be utilized in our case. The Vakhitov- Kolokolov criterion 关19,20兴requiring that the integral N =兰兩u兩2dx be a monotonically growing function of the soliton frequency, was established for stationary solitons. A similar criterion employing Pas a function of Vis valid for solitons moving with a constant velocity 关21–23兴. This latter criterion states that the stability is lost at the point where dP/dV=0, dividing the curve P共V兲into a stable and unstable branch 关8兴. Each point on the curve represents a soliton traveling with some constant velocity and so the points of extremum divide the family of solitons into stable and unstable ones. On the contrary, in our case the oscillatory soliton is visiting the entire P共V兲curve and therefore the analysis of 关8,23兴is not applicable. Finally we show in Sec. VII that the kinetic and canonical soliton momenta are identical, and we analytically calculate the soliton and phonon dispersion curves. II. MODIFIED CONSERVATION LAWS The first three conserved quantities for the unperturbed NLSE are the norm N=兰dx兩u兩2, the momentum P =i 2兰dx共uux ⴱ−uⴱux兲, and the energy E=兰dx关兩ux兩2−兩u兩4兴, where the integration is on −⬁ⱕxⱕ+⬁. The same procedures which are used to prove these conservation laws 关24兴can also be applied to the perturbed NLSE 关Eq. 共1兲兴 in order to obtain the time evolution of N,Pand E: We multiply Eq. 共1兲 by uⴱand its complex conjugate by u, subtract the latter from the former, integrate over xand obtain N ˙=i 冕 −⬁ +⬁ dx共Rⴱu−Ruⴱ兲,共4兲 where the dot denotes the time derivative. Multiplication of Eq. 共1兲by ux ⴱ, addition of the complex conjugate equation multiplied by ux, and integration yields P ˙= 冕 −⬁ +⬁ dx共Rⴱux+Rux ⴱ兲.共5兲 Multiplication of Eq. 共1兲by ut ⴱ, addition of the complex conjugate equation multiplied by ut, and integration yields E ˙=− 冕 −⬁ +⬁ dx共Rⴱut+Rut ⴱ兲.共6兲 Interestingly, for ␤ =0 and time independent force, i.e., R=f共x兲− ␦ u, the r.h.s. of Eq. 共6兲can be written as a time derivative, E ˙=− ⳵ ⳵ t 冕 −⬁ +⬁ dx关fⴱu+fuⴱ− ␦ 兩u兩2兴.共7兲 Thus the perturbed NLSE 关Eq. 共1兲兴 without damping possesses a conserved quantity for arbitrary f共x兲, Etot = 冕 −⬁ +⬁ dx关兩ux兩2−兩u兩4− ␦ 兩u兩2+fⴱu+fuⴱ兴.共8兲 f共x兲can be interpreted as a constant 共external兲force, in contrast to the case of the perturbation R=V共x兲u, where V共x兲 can be understood as a potential 共in which the solitons move兲. In this case the conserved quantity is 关7兴共see also 关25兴兲 Epara tot = 冕 −⬁ +⬁ dx关兩ux兩2−兩u兩4+V共x兲兩u兩2兴.共9兲 In Secs. IV and Vwe will show that a soliton under a constant periodic force f共x兲=aexp共iKx兲performs, on the average, a unidirectional motion, although the average of f vanishes. When we include the damping 共 ␤ ⬎0兲, the time evolution of the total energy is E ˙tot =− ␤ 冕 −⬁ +⬁ dx关2兩ux兩2−4兩u兩4−2 ␦ 兩u兩2+fⴱu+fuⴱ兴. 共10兲 III. COLLECTIVE COORDINATE THEORY For soliton-bearing systems CC theories have been very successful, see Ref. 关26兴for a review and further references. Typically, one starts with the exact one-soliton solution of the unperturbed system and assumes that perturbations of the system change the motion and the shape of the soliton in such a way that both can approximately be described by the one-soliton solution, but with time dependent parameters which are denoted collective coordinates. MERTENS, QUINTERO, AND BISHOP PHYSICAL REVIEW E 81, 016608 共2010兲 016608-2 The 1-soliton solution of the unperturbed NLSE reads 关1兴 u共x,t兲=2i ␩ sech关2 ␩ 共x− ␨ 兲兴e−i共2 ␰ x+ ␾ 兲共11兲 with the real parameters ␩ ⬎0 and ␰ , the soliton position ␨ 共t兲= ␨ 0−4 ␰ t,共12兲 and the phase of the internal oscillation ␾ 共t兲= ␾ 0+4共 ␰ 2− ␩ 2兲t.共13兲 The soliton has amplitude 2 ␩ , width 1/共2 ␩ 兲, velocity V= −4 ␰ and phase velocity of the carrier wave Vph=−2 ␰ +2 ␩ 2/ ␰ . Including the term with ␦ on the rhs of Eq. 共1兲, only the phase is changed: ␾ 共t兲= ␾ 0+关4共 ␰ 2− ␩ 2兲− ␦ 兴t.共14兲 For this reason the term with ␦ need not be treated as a perturbation and will be counted among the unperturbed parts of the NLSE in the following: iut+uxx +2兩u兩2+ ␦ u=f共x,t兲−i ␤ u.共15兲 Therefore we now define the energy as E= 冕 −⬁ +⬁ dx关兩ux兩2−兩u兩4− ␦ 兩u兩2兴,共16兲 which is conserved for the unperturbed NLSE. Using Eq. 共11兲we obtain the soliton energy Esol =16 ␩ ␰ 2−16 3 ␩ 3−4 ␦␩ .共17兲 The definitions of the norm and momentum need not be changed, because all terms with ␦ drop out on the rhs of Eq. 共4兲and 共5兲. The norm and the momentum of the soliton are, respectively, N=4 ␩ ;P=−8 ␩ ␰ .共18兲 Writing P=MV, the soliton mass is related to the norm by M=N/2. A perturbation theory based on the inverse scattering theory 共IST兲results 关1,27兴in a time dependence of ␩ and ␰ and in a change of the simple linear time dependence of ␨ and ␾ in Eqs. 共12兲and 共13兲. Therefore we expect that for sufficiently small perturbations Eq. 共2兲the soliton shape and dynamics can be described by Eq. 共11兲which we take as an ansatz with the four collective coordinates ␩ 共t兲, ␰ 共t兲, ␨ 共t兲and ␾ 共t兲共a similar ansatz has been used for optical solitons where ␩ is related with the beam amplitude, ␰ is a transverse velocity, ␾ is a self-focused beam 关4,28兴兲. We insert Eq. 共11兲into Eqs. 共4兲–共6兲, using f共x,t兲 =aexp关iK共t兲x兴from Eq. 共3兲. Equation 共4兲yields ␩ ˙=−2 ␤ ␩ −a ␲ 2sech Acos B,共19兲 with A共t兲= ␲ 4关K共t兲+2 ␰ 共t兲兴/ ␩ 共t兲,共20兲 B共t兲= ␾ 共t兲+关K共t兲+2 ␰ 共t兲兴 ␨ 共t兲.共21兲 Two of the three terms which result on the rhs of Eq. 共5兲 cancel with the term −8 ␩ ˙ ␰ on the lhs. The remaining term is ␰ ˙=aA sech Acos B.共22兲 Finally, Eq. 共6兲gives 8 ␩ ␰ ˙M−4 ␩ ˙N=0, 共23兲 where Mand Nare some combinations of our collective variables and their derivatives. One solution of Eq. 共23兲 arises by setting both M=0 and N=0; by using an alternative method we will show below that this is the only solution of Eq. 共23兲. Setting M=0 amounts to ␨ ˙=−4 ␰ +a ␲ 2 8 ␩ 2sech Atanh Asin B,共24兲 and letting N=0 is equivalent to ␾ ˙+2 ␨␰ ˙=4共 ␰ 2− ␩ 2兲− ␦ +a ␲ A 2 ␩ sech Atanh Asin B. 共25兲 The four CC equations 共19兲,共22兲,共24兲, and 共25兲can alternatively be obtained by a variational method, because the Lagrangian density for our perturbed NLSE is easily obtained as L=i 2共utuⴱ−ut ⴱu兲−兩ux兩2+兩u兩4+ ␦ 兩u兩2−fuⴱ−fⴱu.共26兲 The dissipative term in the perturbation Eq. 共2兲can be taken into account by generalizing the Euler-Lagrange equation in the following way 关29兴: d dt ⳵ L ⳵ ut ⴱ+d dx ⳵ L ⳵ ux ⴱ− ⳵ L ⳵ uⴱ= ⳵ F ⳵ ut ⴱ,共27兲 with the dissipation function F=i ␤ 共uut ⴱ−uⴱut兲.共28兲 Equation 共27兲is equivalent to our perturbed NLSE. Inserting our CC ansatz and integrating over the system we obtain the CC Lagrangian L=4 ␩ ␾ ˙+8 ␩ ␨␰ ˙−16 ␩ ␰ 2+16 3 ␩ 3+4 ␦␩ −2 ␲ asech Asin B 共29兲 and the CC dissipation function F=− ␤ 共8 ␩ ␾ ˙+16 ␩ ␨␰ ˙兲.共30兲 The CC equations are then obtained by the four generalized Lagrange equations d dt ⳵ L ⳵ ␺ ˙− ⳵ L ⳵ ␺ = ⳵ F ⳵ ␺ ˙,共31兲 where ␺ stands for the four CCs ␩ , ␰ , ␨ , and ␾ . The resulting four ODEs are identical to Eqs. 共19兲,共22兲,共24兲, and 共25兲. NONLINEAR SCHRÖDINGER EQUATION WITH…PHYSICAL REVIEW E 81, 016608 共2010兲 016608-3 This gives an indirect proof that M=0, N=0 provides the only possible solution of Eq. 共23兲. Finally we evaluate Eq. 共8兲, yielding Etot =16 ␩ ␰ 2−16 3 ␩ 3−4 ␦␩ +2 ␲ asech Asin B.共32兲 Here the first three terms are the soliton energy 关Eq. 共17兲兴, while the last term stems from the perturbations. We note that Etot is conserved only in the case of no damping 共 ␤ =0兲and time independent force f共x兲=aexp共iKx兲with constant K, see below Eq. 共7兲. IV. CONSTANT, SPATIALLY PERIODIC FORCE We consider this case because it exhibits some surprising, counterintuitive features: e.g., the soliton position performs oscillations on a length scale that is very different from the spatial period Lof the force. We take a constant Kin Eq. 共3兲, i.e., f共x兲=aexp共iKx兲, and consider only small values of 兩K兩such that the period L =2 ␲ /兩K兩is much larger than the soliton width. We first consider the case with damping 共 ␤ ⬎0兲for which we can expect “steady-state solutions” 关see below Eq. 共36兲兴 for times much larger than a transient time ␶ on the order of 1/ ␤ . The transformation u共x,t兲=⌿共X,t兲exp共iKx兲into a moving frame X=x−Vftwith Vf=2Kleads to the autonomous equation i⌿t+⌿XX +2兩⌿兩2⌿=a+共K2− ␦ 兲⌿−i ␤ ⌿,共33兲 with the nonparametric constant driving term a. Apart from the factor K2− ␦ ªc2, which can be eliminated by scaling time by c2, space by c⬎0, and ⌿by 1/c, Eq. 共33兲is the same as an autonomous equation, which was obtained from the NLSE 关Eq. 共1兲兴 with R= ⑀ exp共i ␻ t兲−i ␤ uby the substitution u=⌿共x,t兲exp共i ␻ t兲, setting ␻ =1 关11,12兴. However, these investigations differ from ours in several respects: In Ref. 关11兴static soliton solutions of Eq. 共33兲were obtained numerically, then an existence and stability chart was constructed on the 共a, ␤ 兲plane. In Ref. 关12兴a singular perturbation expansion was performed at the soliton’s existence threshold. In contrast, we study moving solitons by solving our CC equations, and test the results by simulations, i.e., by numerically solving the NLSE 关Eq. 共1兲兴. Besides the steadystate solutions for ␤ ⬎0共this section兲, for the case ␤ =0 we study stationary solutions and oscillatory solutions, which have a more complicated time dependence 共Sec. V兲. The CC equations from the previous section yield steadystate solutions, in which the driving is compensated by the damping, using the ansatz ␨ = ␨ 0+Vftwith Vf=2Kand constant ␩ , ␰ , and ␾ , uf ⫾共x,t兲=2i ␩ fsech兵2 ␩ f关x− ␨ 共t兲兴其ei共Kx− ␾ f ⫾兲,共34兲 with ␩ f=冑K2− ␦ /2, 共35兲 ␾ f ⫾=⫾ ␲ 2+ arcsin4 ␤ ␩ f ␲ a.共36兲 We denote this as steady-state solutions, because ␾ is constant, in contrast to stationary solutions, where ␾ has a linear time dependence 共Sec. V兲. These solitons have an internal structure due to the factor exp共iKx兲, but no internal oscillations since ␾ f ⫾is constant. In the moving frame these solitons correspond to the above mentioned static solutions of Eq. 共33兲in Ref. 关11兴. We have numerically solved the CC equations for many sets of initial conditions 共IC兲 ␩ 0, ␨ 0, ␰ 0, and ␾ 0. There is a basin of attraction around the solution 共34兲with ␾ f +. The initial conditions always evolve to this solution except when the values of ␩ 0, ␰ 0, and ␾ 0are too far from those of this steady-state solution and when the damping ␤ is too large. An example for this is the parameter set a=0.05, K=0.01, ␦ =−3, ␤ ⱖ ␤ c⬇0.003 with the IC ␩ 0=1, ␰ 0= ␨ 0= ␾ 0=0. Here the soliton vanishes, i.e., its amplitude and energy go to zero while its width goes to infinity. In order to reach the stable steady-state solution one can either reduce ␤ below the critical value ␤ c共which depends on the other parameters and the IC兲, or go closer to ␩ f=0.8660 by reducing ␩ 0by 0.1, for instance, or go closer to ␾ f +by choosing ␾ 0= ␲ /2. In this context it is interesting to consider the total energy 关Eq. 共32兲兴 and its time derivative in which the CC equations can be inserted. We obtain E ˙tot共t兲=− ␤ 兵8 ␩ 关4共 ␰ 2− ␩ 2兲− ␦ 兴+4 ␲ Asech Atanh Asin B其. 共37兲 This is indeed zero for the steady-state solution because ␰ f =−K/2 and thus A⬅0. On the path to that solution, Eq. 共37兲 alternately exhibits both signs. i.e., the total energy performs oscillations which become smaller and smaller while approaching the final value Ef tot. According to Eq. 共35兲, ␦ ⬍K2. As we choose 兩K兩Ⰶ1共see above兲, ␦ is either positive but very small, or ␦ is negative. This results from the CC theory; in Sec. VII we will show that the driven soliton can be stable only for ␦ ⬍0, by taking into account the phonon modes. The quality of the CC theory must be tested by simulations. The soliton shape agrees very well 共Fig. 1兲, but in the simulations the soliton resides on a small constant background because the perturbation f共x兲=aexp共iKx兲does not vanish far away from the soliton. This background is ubg =− a ␻ K eiKx,共38兲 with ␻ K=K2− ␦ −i ␤ ; see also the last term in Eq. 共73兲. Similar to the case of cavity solitons in the DNLSE driven by a biharmonic force 关16兴, here the use of the driving beam f共x兲=aexp共iKx兲supports also cavity solitons sitting on a nonzero background field. Plotting the real and imaginary parts of u, the spatial period L=2 ␲ /兩K兩is observed. The norm density 兩u兩2forms a shelf on which the soliton moves; the shelf height is quantitatively confirmed. The dynamics of the soliton is practically not affected by the background: the time evolution of the soliton position is identical in the CC MERTENS, QUINTERO, AND BISHOP PHYSICAL REVIEW E 81, 016608 共2010兲 016608-4 theory and the simulations 共Fig. 2, right panel兲; only the soliton amplitude differs a little 共left panel兲. All CCs exhibit decreasing oscillations with an intrinsic frequency ⍀which will be discussed in the next section. However, oscillations are not visible in the soliton position ␨ 共t兲in Fig. 2because the linear term dominates the time evolution. By reducing the damping ␤ and by choosing a smaller time scale one can see the oscillations also in ␨ 共t兲, see Fig. 3. In any case, the soliton performs, on the average, a unidirectional motion although the spatial average of the periodic force f共x兲=aexp共iKx兲is zero. Thus, this is a ratchetlike system in which the translational symmetry is broken by the inhomogeneity f共x兲in the NLS equation. The period L=2 ␲ /兩K兩of the inhomogeneity f共x兲is not reflected in the soliton dynamics, because our Eqs. 共1兲–共3兲 could be reduced to the autonomous Eq. 共33兲. For early times 关tⱕ ␶ =O共1/ ␤ 兲兴 the soliton performs the above mentioned oscillations. However, the amplitude of these oscillations is much smaller than L共e.g., L=62.8 for K=0.1兲. When the strength aof the inhomogeneity is strongly increased 共a=1, for instance兲, ␨ 共t兲exhibits a staircase structure. The height of the steps is larger than the amplitude of the above oscillations, but still much smaller than L. However, such values of arepresent a strong perturbation for which the CC theory can no longer be valid. Indeed, in the simulations the soliton soon becomes unstable. V. OSCILLATORY SOLUTIONS In order to study in more detail the intrinsic oscillations in the CCs 共see previous section兲, we consider the case without damping 共 ␤ =0兲. Here oscillatory solutions are possible because the total energy 关Eq. 共32兲兴 is conserved, see Eq. 共37兲. This means that in Eq. 共32兲oscillations of the soliton energy, Eq. 共17兲, are compensated by the oscillations of the term 2 ␲ asech A共t兲sin B共t兲stemming from the perturbations. This is confirmed by inserting numerical solutions of the CC equations into Eq. 共32兲. In order to obtain analytical solutions our approach is to look first for stationary solutions and then to consider small oscillations around them. For the stationary solutions we make the ansatz ␨ = ␨ s+Vst, ␩ = ␩ s, ␰ = ␰ s, and ␾ = ␾ s− ␣ st. The CC equations 共19兲,共22兲,共24兲, and 共25兲yield 0=−a ␲ 2sech Ascos B,共39兲 0=aAssech Ascos B,共40兲 Vs=−4 ␰ s+a ␲ 2 8 ␩ s 2sech Astanh Assin B,共41兲 0 100 200 300 400 500 time 0.7 0.8 0.9 1 1 . 1 η ( t ) 0 100 200 300 400 50 0 time -10 -8 -6 -4 -2 0 ζ(t) (b) ( a) FIG. 2. 共Color online兲The amplitude and position of the soliton obtained from a simulation of the NLSE 共red dashed lines兲and from the numerical solution of the CC equations 共solid lines兲. Parameters: K=−0.01, a=0.05, ␦ =−3, ␤ =0.01, with IC ␰ 0=0, ␨ 0=0, ␾ 0= ␲ /2, and ␩ 0=1. -120 -80 -40 0 x 0 0.4 0.8 1 . 2 |u(x,t*)|2 -200 -100 0 100 20 0 x -0.8 -0.4 0 0.4 Re[u(x,t)]; Im[u(x,t)] (b)(a) FIG. 1. Left panel: soliton moving to the left for tⴱ=250, 500. Simulations of NLSE 共solid lines兲and numerical solutions of CC equations 共dotted lines兲. Right panel: real 共solid line兲and imaginary 共dashed line兲parts of u共x,t兲for t=500. Parameters: K=−0.1, a=0.05, ␦ =−1, ␤ =0.05, with IC ␰ 0=0, ␨ 0=0, ␾ 0=1.69, and ␩ 0=0.5. NONLINEAR SCHRÖDINGER EQUATION WITH…PHYSICAL REVIEW E 81, 016608 共2010兲 016608-5 − ␣ s=4共 ␰ s 2− ␩ s 2兲− ␦ +a ␲ As 2 ␩ s sech Astanh Assin B, 共42兲 with As= ␲ 4共K+2 ␰ s兲/ ␩ s,共43兲 B共t兲= ␾ s+共K+2 ␰ s兲 ␨ s+关共K+2 ␰ s兲Vs− ␣ s兴t.共44兲 Equation 共39兲must be fulfilled for arbitrary t, thus cos B⬅0, which leads to ␣ s=共K+2 ␰ s兲Vs,共45兲 ␾ s+共K+2 ␰ s兲 ␨ s=⫾ ␲ 2,共46兲 sin B=⫾1, 共47兲 respectively. We distinguish two different cases: In case I, ␣ s=0 and we obtain the steady-state solutions of Sec. IV.In the comoving frame these solutions correspond to two exact static solutions of Eq. 共33兲for zero damping 关11,30兴. Here one solution is stable below a critical driving strength, whereas the other one is always unstable. In case II, ␣ s⫽0, which means that we obtain stationary solutions. Here we can restrict ourselves to the case ␨ s= ␰ s =0, because other values yield qualitatively similar results. Using Eqs. 共41兲–共43兲and 共45兲one transcendental equation for ␩ sremains 4 ␩ s 2=− ␦ ⫾ ␲ 2aK 4 ␩ s 2sech Astanh As.共48兲 For the parameter set a=0.05, K=0.1, and ␦ =−3 we get ␩ s =0.866 239 for ␾ s= ␲ /2 and 0.865 811 for ␾ s=− ␲ /2. The numerical solution of the CC equations for these two cases reveals that in the former case we have a stable stationary solution, whereas in the latter case the solution is unstable. For the region around the stable stationary solution we expect that all CCs exhibit oscillations and we assume that these oscillations are harmonic if the amplitudes are sufficiently small. We choose ␾ 0= ␲ /2 and make the ansatz ␨ 共t兲=V ¯ t−a ␨ sin共⍀t兲,共49兲 ␩ 共t兲= ␩ 0+a ␩ 关1 − cos共⍀t兲兴,共50兲 ␰ 共t兲=−a ␰ 关1 − cos共⍀t兲兴,共51兲 ␾ 共t兲= ␾ 0− ␣ t+a ␾ sin共⍀t兲.共52兲 This ansatz takes into account that the soliton first starts its oscillatory motion 共i.e., ␨ −V ¯ tis linear in tfor small t兲, and then the soliton shape changes 共i.e., ␩ is quadratic for small t兲. Since we have considered the stationary solution that belongs to ␰ s= ␨ s=0, we can neglect ␰ compared to Kand obtain for Bin Eq. 共21兲, B= ␾ 0+共a ␾ −Ka ␨ 兲sin共⍀t兲+共KV ¯ − ␣ 兲t.共53兲 Ka ␨ is one order smaller than a ␾ , because 兩K兩Ⰶ1, and we get B= ␾ 0+a ␾ sin共⍀t兲+共KV ¯ − ␣ 兲t.共54兲 Here we can distinguish two limiting cases: in case 1, the linear terms cancel and we have a pure oscillatory behavior. In case 2, the linear term dominates the oscillatory term which can then be neglected. Indeed, taking the second derivative of ␾ with respect to time and using Eqs. 共19兲,共22兲, 共24兲, and 共25兲one obtains ␾ ¨=c1sin共B兲+c2sin共2B兲 +c3cos共B兲+c4cos共2B兲+c5, where ciwith i=1,2,...,5 are functions of the CCs. Considering the above approximations, and in addition when the terms of order of aKV ¯ 共4 ␩ 0 2+ ␦ 兲tare negligible, one realizes that c1,c2,c4, and c5are small compared with c3=⍀2, ⍀2⬇4 ␲ a ␩ 0,共55兲 and hence, the equation for ␾ 共t兲reads ␾ ¨−⍀2cos共 ␾ 兲=0. 共56兲 Integrating this equation once we obtain 0 10203040506070809010 0 time 0.7 0.8 0.9 1 1 . 1 η(t) 0 102030405060708090100 t im e 0 0.5 1 1.5 2 ζ (t) ( b )( a ) FIG. 3. 共Color online兲The amplitude and position of the soliton obtained from a simulation of the NLSE 共red dashed lines兲and from the numerical solution of the CC equations 共solid lines兲. Parameters: K=0.1, a=0.05, ␦ =−3, ␤ =0.001, with IC ␰ 0=0, ␨ 0=0, ␾ 0=0 and ␩ 0 =冑K2− ␦ /2. MERTENS, QUINTERO, AND BISHOP PHYSICAL REVIEW E 81, 016608 共2010兲 016608-6 ␾ ˙=⫾⍀ 冑2冑sin共 ␾ 兲+C,共57兲 where C= ␾ ˙0 2/共2⍀2兲−sin共 ␾ 0兲⬇共共4 ␩ 0 2+ ␦ 兲2兲/共8 ␲ a ␩ 0兲−1. Then, if CⰇ1, the constant term dominates and so ␾ 共t兲goes linearly, i.e., ␾ 共t兲= ␲ /2⫾共4 ␩ 0 2+ ␦ 兲t. Notice that CⰇ1 implies 共4 ␩ 0 2+ ␦ 兲2Ⰷ16 ␲ a ␩ 0, which can be solved for ␩ 0, yielding for a=0.05, ␦ =−3 and ␾ 0= ␲ /2, ␩ 0Ⰶ0.655, or ␩ 0Ⰷ1.077. This is confirmed by numerical solutions of the CC equations which exhibit an oscillatory behavior of ␾ 共t兲 for ␩ 0苸关0.65;1.07兴, and a linear behavior 共plus small oscillations兲outside of this interval. Finally we remark that V ¯ in the oscillatory solutions differs from the velocity Vf=2 Kof the steady-state solutions. A transformation to a frame moving with V⫽2Kwould not simplify the NLSE because terms with ⌿Xwould appear in Eq. 共33兲. VI. STABILITY OF OSCILLATORY SOLUTIONS Our simulations reveal that the driven undamped soliton is stable only for a part of the set of solutions obtained by the CC theory. Naturally we would like to predict, by using the CC theory, which solutions are unstable and to understand what causes the instability. 共The latter point will be discussed in the next section兲. For the NLSE with a general local nonlinearity, but without driving and damping, the Vakhitov-Kolokolov stability criterion 关19,20兴states that solitons are stable if dN/d⌳⬎0. Here Nis the norm and ⌳the so-called spectral parameter in stationary solutions of the form u共x,t兲=⌿共x兲exp共i⌳t兲. However, our oscillatory solutions, Eqs. 共49兲–共52兲inserted into Eq. 共11兲, have a more complicated time dependence than stationary solutions; therefore the criterion cannot be applied here. The same holds for the stability criteria of Barashenkov 关8,23兴which were established for solitons traveling with constant velocity. But here we can get some motivation for how to proceed in our case: Barashenkov showed that dark solitons of the NLSE with generalized nonlinearity are stable if dP ˜ /dV ˜ ⬍0共here the tildes are used to distinguish from our Pand V兲关23兴. This proof was carried over to 共bright and dark兲solitons of the undamped parametrically driven NLSE 关8兴. Here the point dP ˜ /dV ˜ =0 separates a stable from an unstable branch of the curve P ˜ 共V ˜ 兲, but it depends on the type of the solution on which side the stable branch is. For the following it is important to note that this curve P ˜ 共V ˜ 兲represents a family of solutions with different velocities, i.e., each solution is represented by one point on the curve. We make the conjecture that our oscillatory solutions are dynamically unstable, if our curve P共V兲has a branch with negative slope, i.e., dP dV ⬍0, 共58兲 for a finite interval of V. This curve is obtained from its parameter representation P共t兲,V共t兲, where P共t兲=−8 ␩ ␰ is the soliton momentum Eq. 共18兲and the soliton velocity V共t兲= ␨ ˙ is obtained by the rhs of Eq. 共24兲. Each oscillatory solution is represented by its own curve P共V兲. This curve has a finite length, because P共t兲and V共t兲are periodic and remain finite, 关 ␨ 共t兲and ␾ 共t兲, which contain terms linear in t, do not appear in P共t兲, and in V共t兲they appear only via sin B, see Eq. 共24兲兴. Plotting the “stability curve” P共V兲we can immediately see whether there is a branch with negative slope. For the parameter set a=0.05, K=0.1, ␤ =0, and ␦ =−1 and IC ␨ 0= ␰ 0=0 and ␾ 0= ␲ /2, we find a small “stability interval” 0.48ⱕ ␩ 0ⱕ0.52, i.e., an interval of initial conditions for which the solutions are stable. As expected, this interval is situated around the value ␩ s=0.501874 from the stable stationary solution 共48兲. There is another stable regime for ␩ 0ⱖ0.76. When we go far away from the IC for the stationary solution by choosing ␾ 0=0, instead of ␾ 0= ␲ /2, the stability interval around ␩ svanishes. The upper stability regime exists now for ␩ 0ⱖ0.69. When 兩 ␦ 兩is increased, e.g., by choosing ␦ =−3, the stability interval around ␩ s=0.866 239 is much larger than in the case ␦ =−1; for ␾ 0= ␲ /2itis0.7ⱕ ␩ 0ⱕ1.03. The upper stability region is above ␩ 0=1.08. Considering again ␾ 0=0, the stability interval around ␩ sonly shrinks to 0.76ⱕ ␩ 0ⱕ0.97, but does not vanish because it was much larger than in the case ␦ =−1. The upper stable region is above ␩ 0=1.02. Thus the conclusion is that an increase of 兩 ␦ 兩widens the regions with stable soliton solutions. All the stability regions given above are confirmed by simulations for the perturbed NLSE, with an error of less than 1%. At some of the boundaries of the above stability intervals there is a drastic change in the shape of the solutions of the CC equations and the stability curve. E.g., if we choose ␩ 0 =0.75 共Fig. 4兲, which is just below the stability regime ␩ 0 ⱖ0.76 共see above兲, we obtain very anharmonic oscillations in all CCs and the stability curve has a long branch with negative slope 共Fig. 5兲. The simulations indeed show that the soliton becomes unstable very quickly and vanishes 共Fig. 4兲. However, a slight change of ␩ 0to the value 0.76 produces a stability curve that has only one branch with a positive slope 共Fig. 5兲. Here the soliton is indeed stable in the simulations 共Fig. 6兲. At the boundaries of other stability intervals the change of the shape of the solutions is less dramatic. E.g., for the case ␦ =−1 and ␩ 0=0.46 the oscillations in the CCs are small and harmonic. The stability curve has only a very short branch with negative slope which is visited in the time evolution only for very short time intervals 共Fig. 7兲. Here the soliton is indeed stable for a relatively long time 共Fig. 8兲. This soliton lifetime is increasingly reduced when ␩ 0is reduced by which the negative-slope branch in P共V兲becomes longer. Thus this curve predicts not only whether the soliton is stable, but also gives an estimate for its lifetime when it is unstable. VII. SOLITON AND PHONON DISPERSION CURVES For the soliton dispersion curve Esol共P兲we need the canonical momentum Pof the soliton. Pmust fulfill the Hamilton equation ␨ ˙= ⳵ H ⳵ P,共59兲 where ␨ is the soliton position. We start from the CC Lagrangian Eq. 共29兲and obtain NONLINEAR SCHRÖDINGER EQUATION WITH…PHYSICAL REVIEW E 81, 016608 共2010兲 016608-7 ⳵ L ⳵ ␾ ˙=4 ␩ =N,共60兲 as the angular momentum of the internal oscillation of the soliton, and ⳵ L ⳵␰ ˙=8 ␩ ␨ =2N ␨ ªX,共61兲 as the variable which is canonically conjugate to ␰ . Using a Legendre transformation we obtain the Hamilton function H=X ␰ ˙+N ␾ ˙−L,共62兲 H共X, ␰ ;N, ␾ 兲=4N ␰ 2−1 12N3− ␦ N+2 ␲ asech Asin B, 共63兲 with A=2 ␲ K+2 ␰ 2N,B= ␾ +K+2 ␰ 2NX.共64兲 We now make an ansatz for a canonical transformation to the following new set of variables: P=−2N ␰ , ␨ =1 2NX,共65兲 0 20406080100 t im e 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 . 8 η ( t ) 0 2040608010 0 t im e -35 -30 -25 -20 -15 -10 -5 0 5 ζ(t) 0 2040608010 0 time -0.05 0 0.05 0.1 0.15 ξ (t) 0 2040608010 0 t im e -2 0 2 4 6 8 φ(t) (b)(a) (c) (d) FIG. 4. The evolution of ␩ , ␨ , ␰ , and ␾ obtained from the numerical solutions of the CC equations 共solid lines兲and the evolution of ␩ and ␨ from a simulation of the NLSE 共dashed lines兲. Parameters: K=0.1, a=0.05, ␦ =−1, ␤ =0, with IC ␰ 0=0, ␨ 0=0, ␾ 0= ␲ /2, and ␩ 0 =0.75. -0.4 -0.3 -0.2 -0.1 0 0.1 0. 2 v -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 P -0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0 0.0 2 v -0.08 -0.06 -0.04 -0.02 0 0.02 P (b)(a) FIG. 5. “Stability curve” Pversus V. Left panel: ␩ 0=0.75. Right panel ␩ 0=0.76. Other parameters as in Fig. 4. MERTENS, QUINTERO, AND BISHOP PHYSICAL REVIEW E 81, 016608 共2010兲 016608-8 N ˜ =N, ␾ ˜ = ␾ +g共N, ␰ ,X兲.共66兲 This means that Pis identical to the kinetic momentum in Eq. 共18兲,N ˜ is chosen to be equal to Nand gstill has to be determined such that the four fundamental Poisson brackets are fulfilled. This yields g= ␰ X/N. The new Hamiltonian reads H共P, ␨ ;N, ␾ ˜ 兲=1 NP2−1 12N3− ␦ N+2 ␲ asech Asin B, 共67兲 with A= ␲ N 冉 K−P N 冊 ,B= ␾ ˜ +K ␨ .共68兲 The first two terms in Eq. 共67兲agree with literature results on the unperturbed NLSE 关3兴. We have checked that the four Hamiltonian equations which result from Eq. 共67兲are indeed equivalent to the four CC equations in Sec. III for the case without damping. The above results hold for the force f共x,t兲 in Eq. 共3兲with arbitrary K共t兲. In the following we return to the case of constant Kin order to calculate analytically the stability and dispersion curves wherever it is possible. Let us consider the stability intervals around the stationary solutions as given in Sec. VI. Except for the close vicinity to the boundaries of the stability intervals, the oscillations in the CCs are nearly harmonic and can be well approximated by the Eqs. 共49兲–共52兲. Inserting into P=−8 ␩ ␰ , neglecting the second harmonic, and using V= ␨ ˙, one can easily see that P共V兲is a straight line with slope dP dV =8 ␩ 0a ␰ ⍀a ␨ ⬎0, 共69兲 because for ␩ 0⬍ ␩ sboth a ␰ and a ␨ are positive and for ␩ 0⬎ ␩ sboth are negative. The latter also holds for the upper stability intervals ␩ 0ⱖ0.76 for ␦ =−1 and ␩ 0ⱖ1.08 for ␦ =−3. Near or at the boundaries of the stability intervals, the oscillations in the CCs are very anharmonic and therefore the calculation leading to Eq. 共69兲is not possible. P共V兲is a curved line which can be calculated using the numerical solution of the CC equations. We now turn to the soliton dispersion curve Esol共P兲. Both P=−8 ␩ ␰ and Esol in Eq. 共17兲consist of powers of ␩ and ␰ . 0 2040608010 0 t im e 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 η(t) 0 2040608010 0 t im e -8 -6 -4 -2 0 ζ (t) 0 2040608010 0 t im e -0.005 0 0.005 0.01 0.015 0.02 0.025 ξ (t) 0 2040608010 0 t im e -60 -50 -40 -30 -20 -10 0 φ(t) (b) ( a) ( c) (d) FIG. 6. The evolution of ␩ , ␨ , ␰ , and ␾ obtained from the numerical solutions of the CC equations 共solid lines兲and the evolution of ␩ and ␨ from a simulation of the NLSE 共dashed lines兲. Parameters: K=0.1, a=0.05, ␦ =−1, ␤ =0, with IC ␰ 0=0, ␨ 0=0, ␾ 0= ␲ /2, and ␩ 0 =0.76. The simulations were carried forward to a final time tf=1000. 0.048 0.05 0.052 0.054 0.05 6 v 0 0.01 0.02 0.03 0 . 04 P FIG. 7. “Stability curve” close to the boundary of an unstable regime. Parameters: K=0.1, a=0.05, ␦ =−1, ␤ =0, with IC ␰ 0=0, ␨ 0=0, ␾ 0= ␲ /2, and ␩ 0=0.46. NONLINEAR SCHRÖDINGER EQUATION WITH…PHYSICAL REVIEW E 81, 016608 共2010兲 016608-9