Forced nonlinear Schrödinger equation with arbitrary nonlinearity

PHYSICAL REVIEW E 85, 046607 (2012) Forced nonlinear Schr¨ odinger equation with arbitrary nonlinearity Fred Cooper,1,2,*Avinash Khare,3,†Niurka R. Quintero,4,‡Franz G. Mertens,5,§and Avadh Saxena2, 1Santa Fe Institute, Santa Fe, New Mexico 87501, USA 2Theoretical Division and Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA 3Indian Institute of Science Education and Research, Pune 411021, India 4IMUS and Departamento de Fisica Aplicada I, EUP Universidad de Sevilla, 41011 Sevilla, Spain 5Physikalisches Institut, Universit¨ at Bayreuth, D-95440 Bayreuth, Germany (Received 29 November 2011; published 24 April 2012) We consider the nonlinear Schr¨ odinger equation (NLSE) in 1 +1 dimension with scalar-scalar self-interaction g2 κ+1(ψψ)κ+1in the presence of the external forcing terms of the form re−i(kx+θ)−δψ. We find new exact solutions for this problem and show that the solitary wave momentum is conserved in a moving frame where vk=2k. These new exact solutions reduce to the constant phase solutions of the unforced problem when r→0. In particular we study the behavior of solitary wave solutions in the presence of these external forces in a variational approximation which allows the position, momentum, width, and phase of these waves to vary in time. We show that the stationary solutions of the variational equations include a solution close to the exact one and we study small oscillations around all the stationary solutions. We postulate that the dynamical condition for instability is that dp(t)/d ˙ q(t)<0, where p(t) is the normalized canonical momentum p(t)=1 M(t) ∂L ∂˙ q,and ˙ q(t) is the solitary wave velocity. Here M(t)=dxψ(x,t)ψ(x,t). Stability is also studied using a “phase portrait” of the soliton, where its dynamics is represented by two-dimensional projections of its trajectory in the four-dimensional space of collective coordinates. The criterion for stability of a soliton is that its trajectory is a closed single curve with a positive sense of rotation around a fixed point. We investigate the accuracy of our variational approximation and these criteria using numerical simulations of the NLSE. We find that our criteria work quite well when the magnitude of the forcing term is small compared to the amplitude of the unforced solitary wave. In this regime the variational approximation captures quite well the behavior of the solitary wave. DOI: 10.1103/PhysRevE.85.046607 PACS number(s): 05.45.Yv, 11.10.Lm, 63.20.Pw I. INTRODUCTION The nonlinear Schr¨ odinger equation (NLSE), with cubic and higher nonlinearity counterparts, is ubiquitous in a variety of physical contexts. It has found several applications, including in nonlinear optics where it describes pulse propagation in double-doped optical fibers [1] and in Bragg gratings [2], and in Bose-Einstein condensates (BECs) where it models condensates with two- and three-body interactions [3,4]. Higher order nonlinearities are found in the context of Bose gases with hard core interactions [5] and low-dimensional BECs in which quintic nonlinearities model three-body interactions [6]. In nonlinear optics a cubic-quintic NLSE is used as a model for photonic crystals [7]. Therefore it is important to ask the question how will the behavior of these systems change if a forcing term is also included in the NLSE. The forced nonlinear Schr¨ odinger equation (FNLSE) for an interaction of the form (ψψ)2has been recently studied [8,9] using collective coordinate (CC) methods such as timedependent variational methods and the generalized traveling wave method (GTWM) [10]. In Refs. [8,9] approximate stationary solutions to the variational solution were found and a criterion for the stability of these solutions under small perturbations was developed and compared to numerical *[email protected] †[email protected] ‡[email protected] §[email protected] a[email protected]v simulations of the FNLSE. Here we will generalize our previous study to arbitrary nonlinearity (ψψ)κ+1, with a special emphasis on the case κ=1/2. That is, the form of the FNLSE we will consider is i∂ ∂tψ+∂2 ∂x2ψ+g(ψψ)κψ+δψ =re−i(kx+θ).(1.1) The parameter rcorresponds to a plane wave driving term. The parameter δarises in discrete versions of the NLSE used to model discrete solitons in optical wave guide arrays and is a cavity detuning parameter [11]. We will find that having δ<0 allows for constant phase solutions of the CC equations. The externally driven NLSE arises in many physical situations such as charge density waves [12], long Josephson junctions [13], optical fibers [14], and plasmas driven by rf fields [15]. What we would like to demonstrate here is that the stability criterion for the FNLSE solitons found for κ=1 works for arbitrary κ⩽2, and that the collective coordinate method works well in predicting the behavior of the solitary waves when the forcing parameter ris small compared to the amplitude of the unforced solitary wave. The paper is organized as follows. In Sec. II we show that in a comoving frame where y=x+2kt, the total momentum of the solitary wave Pvas well as the energy of the solitary wave Evis conserved. In Sec. III we review the exact solitary waves for r=0 and κarbitrary. We show using Derrick’s theorem [16] that these solutions are unstable for κ>2 and arbitrary δ, which we later verify in our numerical simulations. In Sec. IV we find exact solutions to the forced problem for r= 0 and find both plane wave as well as solitary wave and 046607-1 1539-3755/2012/85(4)/046607(24) ©2012 American Physical Society COOPER, KHARE, QUINTERO, MERTENS, AND SAXENA PHYSICAL REVIEW E 85, 046607 (2012) periodic solutions for arbitrary κ. We focus mostly on the case κ=1/2. We find both finite energy density as well as finite energy solutions. In Sec. V we discuss the collective coordinate approach in the laboratory frame. We will use the form of the exact solution for the unforced problem with time dependent coefficients as the variational ansatz for the traveling wave for κ<2. This is a particular example of a collective coordinate (CC) approach [8,17–19]. We will assume the forcing term is of the form re−i(kx+θ)−δψ. We will choose four collective coordinates, the width parameter β(t), the position q(t) and momentum p(t), as well as the phase φ(t)of the solitary wave. These CCs are related by the conservation of momentum in the comoving frame. We will derive the effective Lagrangian for the collective coordinates and determine the equations of motion for arbitrary nonlinearity parameter κ.In Sec. VI we show that the equations of motion for the collective coordinates simplify in the comoving frame. For arbitrary κwe determine the equations of motion for the collective coordinates, the stationary solutions (β=˙ q=φ=const) as well as the linear stability of these stationary solutions. We then specialize to the case κ=1/2, and discuss the linear stability of the stationary solutions. The real or complex solutions to the small oscillation problem give a local indication of stability of these stationary solutions. This analysis will be confirmed in our numerical simulations. For κ=1/2 and rA, where A is the amplitude of the solitary wave, we find in general three stationary solutions. Two are near the solutions for r=0, one being stable and another unstable, and one is of much smaller amplitude but turns out to be a stable solitary wave. A more general question of stability for initial conditions having an arbitrary value of β(t=0) ≡β0is provided by the phase portrait of the system found by plotting the trajectories of the imaginary vs the real part of the variational wave function starting with an initial value of β0. These trajectories are closed orbits as shown in Ref. [9]. Stability is related to whether the orbits show a positive (stable) or negative (unstable) sense of rotation or a mixture (unstable). Another method for discussing stability is to use the dynamical criterion used previously [9] in the study of the case κ=1, namely whether the p(˙ q) curve has a branch with negative slope. If this is true, this implies instability. Here p(t) is the normalized canonical momentum p(t)=1 M(t) ∂L ∂˙ q,M(t)=dxψ(x,t)ψ(x,t) and ˙ q(t)=v(t)is the velocity of the solitary wave. These two approaches (phase portrait and the slope of the p(v) curve) give complementary approaches to understanding the behavior of the numerical solutions of the FNLSE. In Sec. VII we discuss how damping modifies the equations of motion by including a dissipation function. We find that damping only effects the equation of motion for βamong the CC equations. This damping allows the numerical simulations to find stable solitary wave solutions. In Sec. VIII we discuss our methodology for solving the FNLSE. We explain how we extract the parameters associated with the collective coordinates from our simulations. We first show that our simulations reproduce known results for the unforced problem as well as show that the exact solutions of the forced problem we found are metastable. On the other hand the linearly stable stationary solutions to the CC equations are close to exact numerical solitary waves of the forced NLSE with time independent widths, only showing small oscillations about a constant value of β. In the numerical solutions of both the PDEs and the CC equations we establish two results. First, for small values of the forcing term rthe CC equations give an accurate representation of the behavior of the width, position, and phase of the solitary wave determined by numerically solving the NLSE. Second, both the phase portrait and p(v) curves allow us to predict the stability or instability of solitary waves that start initially as an approximate variational solitary wave of the form of the exact solution to the unforced problem. In Sec. IX we state our main conclusions. II. FORCED NONLINEAR SCHR ¨ ODINGER EQUATION (FNLSE) The action for the FNLSE is given by Ldt =dtdxiψ∂ ∂tψ−ψ xψx−g κ+1(ψψ)κ+1 −δψψ+rei(kx+θ)ψ+re−i(kx+θ)ψ.(2.1) As shown in Ref. [8] the energy E=dxψ xψx−g κ+1(ψψ)κ+1−δψψ +rei(kx+θ)ψ+re−i(kx+θ)ψ(2.2) is conserved. Varying the action, Eq. (2.1) leads to the equation of motion: i∂ ∂tψ+∂2 ∂x2ψ+g(ψψ)κψ+δψ =re−i(kx+θ).(2.3) We notice that the equation of motion is invariant under the joint transformation r→−r,ψ→−ψ, thus if ψ(x,t,r)isa solution, then so is −ψ(x,t,−r) a solution. Letting ψ(x,t)= e−i(kx+θ)u(x,t) we obtain the autonomous equation i∂u(x,t) ∂t −2k∂u(x,t) ∂x −(k2−δ)u(x,t) +∂2u(x,t) ∂x2+g(uu)κu=r. (2.4) Changing variables from xto ywhere y=x+2kt =x+vkt we have for u(x,t)→v(y,t) i∂v(y,t) ∂t −(k2−δ)v(y,t)+∂2v(y,t) ∂y2+g(vv)κv=r. (2.5) Note that with our conventions, the mass in the Schr¨ odinger equation obeys 2m=1, or m=1/2 so that k=mvk=vk/2. This equation in the moving frame can be derived from a related action Sv=dtdyiv∂ ∂tv−v yvy−g κ+1(vv)κ+1 +(k2−δ)vv+rv +rv.(2.6) Multiplying Eq. (2.5) on the left by v yand adding the complex conjugate, we get an equation for the time evolution of the momentum density in the moving frame: ρv=i 2(vv y−vvy), 046607-2 FORCED NONLINEAR SCHR ¨ ODINGER EQUATION WITH ... PHYSICAL REVIEW E 85, 046607 (2012) namely ∂ρv ∂t +∂j(y,t) ∂y =rv y+rvy.(2.7) Here j(y,t)=i 2(vvt−vv t)+|vy|2 +(δ−k2)|v|2+g κ+1|v|2κ+2.(2.8) Integrating over all space we get d dt dyρv(y,t)=F[y=+∞]−F[y=−∞],(2.9) where F[y]=−j(y,t)+rv(y,t)+rv(y,t).(2.10) If the value of the boundary term is the same (or zero) at y=±∞then we find that the momentum Pvin the moving frame is conserved: Pv=const =dy i 2(vv y−vvy).(2.11) Using the fact that ψ=ue−i(kx+θ)and defining P=dx i 2(ψψ x−ψψx),(2.12) we obtain the conservation law P(t)+M(t)k=Pv=const,(2.13) where M(t)=dxψψ. If we further define p(t)=P(t) M(t),(2.14) then we get the relationship Pv=M(t)[p(t)+k]=const.(2.15) This equation will be useful when we consider variational approximations for the solution. The second action Eq. (2.6) also leads to a conserved energy: Ev=dxv yvy−g κ+1(vv)κ+1 +(k2−δ)vv+rv +rv.(2.16) Using the connection that ψ(x,t)=e−i(kx+θ)v(y,t) we find that Eand Evare related. Using the fact that ψ∗ xψx=k2v∗v+v∗ yvy+ik(v∗vy−vv∗ y),(2.17) we obtain E=Ev−2kPv.(2.18) III. EXACT SOLITARY WAVE SOLUTIONS WHEN r=0 Before discussing the exact solutions to the forced NLSE, let us review the exact solutions for r=0 since these will be used as our variational trial functions later in the paper. We can obtain the exact solitary wave solutions when r=0as follows. We let ψ(x,t)=Asechγ[β(x−vt)]ei[p(x−vt)+φ(t)].(3.1) Demanding that we have a solution by matching powers of sech we find p=v/2,γ=1/κ, A2κ=β2(κ+1) gκ2, (3.2) φ=(p2+β2/κ2+δ)t+φ0. It is useful to connect the amplitude Ato βand the mass Mof the solitary wave. M=A2dxψ∗(βx)ψ(βx)=A2 βC0 1; (3.3) C0 1=+∞ −∞ dy sech2/k(y)=√π1 κ 1 2+1 κ. One can show that the energy E=+∞ −∞ dx(x) of the solitary waves is given by E=√π1 κ κ+2 2β2−κ(κ+1) gκ21/κ p2−δ+β2(κ−2) κ2(κ+2) . (3.4) For δ=0, κ>2 these solitary waves are unstable [18,19]. For κ=2 the solitary wave is a critical one. In this case, the energy is independent of the width of the solitary wave. Further, in the rest frame v=0, the energy is zero when δ=0. When δis not zero there is a special constant phase solution of the form Eq. (3.1). The condition for φto be independent of tis p2+β2/κ2+δ=0,(3.5) which is only possible for negative δ. For this solution, when p=0 we find the relationship β2=−κ2δ. (3.6) This solution will be the r=0 limit of some of the solutions we will find below. A. Stability when r=0 For the unforced NLSE we can use the scaling argument of Derrick [16] to determine if the solutions are unstable to scale transformation. We have discussed this argument in our paper on the nonlinear Dirac equation [20]. The Hamiltonian is given by H=dx 1 2mψ xψx−δψψ−g κ+1(ψψ)κ+1.(3.7) It is well known that using stability with respect to scale transformation to understand domains of stability applies to this type of Hamiltonian. This Hamiltonian can be written as H=H1−δH2−H3,(3.8) where Hi>0. Here δcan have either sign. If we make a scale transformation on the solution which preserves the mass M=ψψdx, ψβ→β1/2ψ(βx),(3.9) 046607-3 COOPER, KHARE, QUINTERO, MERTENS, AND SAXENA PHYSICAL REVIEW E 85, 046607 (2012) we obtain H=β2H1−δH2−βκH3.(3.10) The first derivative is ∂H/∂β =2βH1−κβκ−1H3. Setting the derivative to zero at β=1 gives an equation consistent with the equations of motion: κH3=2H1.(3.11) The second derivative at β=1 can now be written as ∂2H ∂β2=κ(2 −κ)H3.(3.12) The solution is therefore unstable to scale transformations when κ>2. This result is independent of δ. However, once one adds forcing terms, it is known from the study of the κ=1 case [8] that the windows of stability as determined by the stability curve p(v) as well as by simulation of the FNLSE increase as δis chosen to be more negative. In those simulations the two methods agreed to within 1%. IV. EXACT SOLUTIONS OF THE FORCED NLSE FOR r= 0 For r=0 one can have time dependent phases in the traveling wave solutions, however when r= 0 one is restricted to looking for traveling wave solutions with time independent phases. That is, if we consider a solution of the form ψ(x,t)=e−i(kx+θ)f(y); y=x+2kt, (4.1) then fsatisfies f(y)−k2f+g(ff)κf=r;k2=k2−δ. (4.2) A. Plane wave solutions First let us consider the plane wave solution of Eq. (2.3): ψ(x,t)=aexp [−i(kx +θ)],(4.3) or equivalently the constant solution f=ato Eq. (4.2).This is a solution provided asatisfies the equation (here acan be positive or negative) r=g(a2)κa−(k2−δ)a. (4.4) This solution has finite energy density, but not finite energy in general. The energy density (x) is given by (x)=g(2κ+1)a2(κ+1) κ+1−(k2−δ)a2.(4.5) Since Pv=0, the energy is the same in the comoving frame. This is the lowest energy solution for unrestricted g,a,k2, and δ. The solitary wave solutions discussed below will have finite energy with respect to this “ground state” energy. There is a special zero energy solution that is important. This solution has the restriction that g(2κ+1)a2κ κ+1=(k2−δ).(4.6) B. 2κ=integer When 2κ=integer then one also can simply find solutions. The differential equation is f −(k2)f+g(f∗f)κf=r. (4.7) Note that again this equation is invariant under the combined transformation f→−fand r→−r. Thus we look for solutions of the form f±=±[a+u(y)]with a>0; u(y)>0 so that for this ansatz we have |a+u(y)|=a+u(y).(4.8) It is sufficient to consider f+=a+u(y) and generate the second solution by symmetry (r→−r;a→−a;u→−u). For f+we have letting 2κ+1=N=integer r=−k2a+gaN,(4.9) and u =(k2−gNaN−1)u−g N  m=2 umaN−mN m.(4.10) Integrating once and setting the integration constant to zero to keep the energy of the solitary wave relative to a plane wave finite, we obtain (u)2=(k2−gNaN−1)u2−g N  m=2 2um+1aN−mN! (m+1)!(N−m)!. (4.11) This equation will have a solution as long as α=k2− gNaN−1>0. We now discuss solutions to this equation when N=2(κ=1/2) and N=3(κ=1). C. κ=1/2 For κ=1/2, N=2, the equation we need to solve for f+ is (u)2=αu2−2gau3−gu4/2; α=k2−2ga. (4.12) This equation has a solution of the form u(y)=bsech1/κ (βy),(4.13) so that f(y) has a solution of the form f(y)=a+bsech2β(y,t),(4.14) where both a,b > 0,and fobeys the equation f −(k2)f+g|f|f=r. (4.15) The wave function is given by ψ(x,t)=exp[−i(kx +θ)][a+bsech2β(x+2kt)].(4.16) First consider the case that a>0 and b>0, so that |f|=f. Substituting the ansatz Eq. (4.14) into Eq. (4.15) we obtain a2g+2abg sech2(βy)−a(k)2+b2gsech4(βy) −b(k)2sech2(βy)−6bβ2sech4(βy) +4bβ2sech2(βy)−r=0.(4.17) 046607-4 FORCED NONLINEAR SCHR ¨ ODINGER EQUATION WITH ... PHYSICAL REVIEW E 85, 046607 (2012) Equating powers of sech we obtain the conditions −k2a+ga2=r;4β2−k2+2ga =0; (4.18) −6bβ2+gb2=0, or, equivalently, a=k2−k4+4rg 2g,β 2=1 4k4+4rg, b =6β2 g. (4.19) The assumption a>0 and b>0 requires for consistency that k2>0 and r<0. So we can rewrite a=k2−k4−4|r|g 2g, (4.20) β2=1 4k4−4|r|g, b =6β2 g. We also have the restriction that α=k2−2ga > 0.(4.21) The energy density corresponding to this solution is easily calculated and is given by (x)=4ga3 3−(k2−δ)a2+b2[k2−δ−2ga +(2/3)bg] ×sech4[β(x+2kt)] −(4/3)gb3sech6[β(x+2kt)]. (4.22) Observe that the constant term is exactly the same as the energy density of the solution (4.3) at κ=1/2. Hence the energy of the pulse solution [over and above that of the solution (4.3)]is given by E=384β5 5g2.(4.23) Because of symmetry there is also a solution with f→ −f,r →−rso that we can write the two solutions connected by this symmetry as follows: f(y)=−sgn(r)[a+bsech2(βy)] (4.24) with a,β, and bgiven by Eq. (4.20). Thus the wave function can be written as ψ(x,t)=−sgn(r)[a+bsech2β(x+2kt)]e−i(kx+θ).(4.25) Notice that the boundary conditions (BCs) on this solution are that for ψ(x)at±∞ the solution goes to the plane wave solution ψ(x=±∞)→−sgn(r)aexp[−i(kx +θ)].(4.26) Thus the BC on solving this equation numerically is the mixed boundary condition ikψ(x=±∞,t)+ψx(x=±∞,t)=0.(4.27) On the other hand, if we go into the yframe where u(y=±∞,t)→−sgn(r)a, (4.28) then the BC for the uequation is uy(y=±∞,t)=0.(4.29) FIG. 1. (Color online) Exact f(y) (red dashed line) versus variational solution f2+(y) (black solid line) for k=−0.1, δ=−1, g=2, r=−0.01. Consider the case where g=2,k =−0.1,δ =−1, r=−0.01. For this choice of parameters the exact solution is given by f(y)=0.727 191 sech2(0.492 338y)+0.010 103 1.(4.30) The stationary solution found using our variational ansatz which we will discuss below [see Eq. (6.65)] yields the approximate solution: f2+(y)=0.745 042 sech2(0.498 345y),(4.31) which is a reasonable representation of the exact solution as seen in Fig. 1. We will show later by numerical simulation, that this solution is metastable in that βremains constant for a short period of time and then oscillates. In the numerical simulation of Fig. 9the parameters used in the simulation are κ=1/2,δ =−1,g =2,k =0.1, r=−0.075. For that case, f(y)=0.486 113 sech2(0.402 539y)+0.090 462 2.(4.32) Note that here ris 14% of Aso that the variational solution is worse. The unstable stationary variational solution corresponding to this is f2+(y)=0.745 535 sech2(0.498 509y).(4.33) To make a comparison we plot the square of these solutions in Fig. 2. FIG. 2. (Color online) Exact f2(y) (red dashed line) versus variational solution f2 2+(y) (black solid line) for δ=−1, g=2, k=0.1, r=−0.075. These parameters were used in the simulations shown in Fig. 9. 046607-5 COOPER, KHARE, QUINTERO, MERTENS, AND SAXENA PHYSICAL REVIEW E 85, 046607 (2012) FIG. 3. Finite energy solution for κ=1/2, β=1, normalized to unit height. D. Finite energy solutions From Eq. (4.22) we notice that we can have a finite energy solution [with energy as given by Eq. (4.23)] when the constant term in the energy density is zero. This leads to the condition 4ga =3k2, and yields b=−a. This solution is not asmall perturbation on the r=0 solution which has a=0,b = 0. Instead for the finite energy solution a=−band the form of the solution is f(y)=Atanh2βy, (4.34) which has the appearance of Fig. 3. If we insert the solution Eqs. (4.34) into (4.15), we obtain Ag|A|tanh4(βy)+2Aβ2−Ak2tanh2(βy) +6Aβ2tanh4(βy)−8Aβ2tanh2(βy)−r=0.(4.35) Thus Eq. (4.34) is a solution provided r=2Aβ2;β2=−k2/8; |A|=−6β2/g, (4.36) which requires that both g<0 and k2<0. We see that this solution has the symmetry that rchanges sign with Aand thus the sign of fdepends on the sign of r. Because the coupling constant needs to be negative, the quantum version of this theory would not have a stable vacuum. We can write this solution as f= r |r| 6β2 |g|tanh2βy, (4.37) where now β2=|k2|/8 and rhas the special values r±= ±(3/16)(k2)2/|g|. 1. Periodic solutions for κ=1/2 For κ=1/2 one can easily verify that there are periodic solutions of Eq. (4.2) of the form f+(y)=bdn 2(βy,m)+a, (4.38) where dn(x,m) is the Jacobi elliptic function (JEF) with modulus m.Againfora>0,b > 0wehave|f+|=f+. Matching coefficients of powers of dn one obtains a=√1−m+m2k2−(2 −m)4gr +k4 2g√1−m+m2, (4.39) b=6β2 g,β 2=k4+4gr 4√1−m+m2. The requirement that a>0 translates into a restriction on r, namely r<−3(1 −m)k4/4g(2 −m)2. E. κ=1 When κ=1, the equation for fis f −k2f+g|f|2f=r. (4.40) Again it is only necessary to consider f+=[a+u(y)] with both a>0 and u(y)>0. For κ=1, N=3, the equation we need to solve for u(y)is (u)2=αu2−2gau3−gu4/2.(4.41) This equation has two solutions of the form u(y)=b c±cosh(βy)(4.42) (note that y=x+2kt). These solutions were found earlier by Barashenkov and collaborators [21,22]. Our condition on ufor f+requires that b>0,c > 0 and we choose the +solution. This leads to ψ+(x,t)=exp[−i(kx +θ)]a+b c+cosh β(x+2kt), (4.43) provided β=√α;b=√2α 2g2a2+gα; (4.44) c=√2ag 2a2g2+gα;α=k2−3ga2. Since α>0, we require k2>3ga2>0. Here r=−k2a+ ga3<−2ga3is negative. When we let a→−a,b →−b, c→−c, then rchanges sign and becomes positive. Thus we can write ψ±(x,t)=−sgn(r)exp[−i(kx +θ)] ×|a|+ |b| |c|+cosh β(x+2kt).(4.45) These solutions are nonsingular since 1−c2=gα/[gα +2g2a2]>0.(4.46) The energy density corresponding to these solutions is given by (y)=3ga4 2−k2a2+2b2(k2−3ga2) [c+cosh(βy)]2−2b2(cβ2+gab) [c+cosh(βy)]3 +b2[β2(c2−1) −(g/2)b2] [c+cosh(βy)]4.(4.47) Note that the constant term is exactly the same as the energy of the solution (4.3) as given by Eq. (4.5) at κ=1. Hence the energy of the κ=1 pulse solution [over and above that of the solution (4.3)] is again finite and given by E=8α5/2 gα +2a2g2I2−16√2α5/2ag [gα +2a2g2]3/2I3−8gα7/2 [gα +2a2g2]2I4, (4.48) 046607-6 FORCED NONLINEAR SCHR ¨ ODINGER EQUATION WITH ... PHYSICAL REVIEW E 85, 046607 (2012) where Ij=∞ 0 dy [c+cosh(βy)]j,j=2,3,4.(4.49) We have I1=∞ 0 dx a+bcosh(x)=2 √b2−a2tan−1√b2−a2 b+a, if b2>a 2.(4.50) Note that in our case a=c,b =1,c2<1. From here it is easy to calculate the integrals I2,I3,I4and hence show that the energy of the +and the −solution [over and above the solution (4.3)]isgivenby E=8α1/2 3g(α+3a2g)−8√2a √g(α+2a2g) ×tan−11+2a2g α−2a2g α.(4.51) 1. Periodic solutions for κ=1 For κ=1 we can generalize the solitary wave solution for f+, f −k2f+gf 3=r, (4.52) namely f=a+b c+cosh[dy]=(ac +b) sech[dy]+a 1+csech[dy],(4.53) and obtain a periodic solution in terms of the Jacobi elliptic functions. The generalization of sech[dy] is the Jacobi elliptic function dn(dy,m) which also has the property |dn|= dn(dy,m). One finds that f+=a+hdn(dy,m) 1+cdn(dy,m),(4.54) with a>0,h > 0,c > 0 obeys |f+|=f+and is an exact solution to Eq. (4.52) provided that r=(ac +h)(agh −ck2) 2c2;d2=−3agh −ck2 2c(m−2) ,(4.55) h2=4d2 g[2 +αga2(1 −k2α)], (4.56) (1 −m)h2=1+2αga2−k2α 4d2α3g2a2. Here α=c/gah which obeys a cubic equation, 16(1 −m)α3d4ga2 =[2 +αga2(1 −k2α)][1 +2αga2−k2α].(4.57) F. κ=3/2 For κ=3/2,N =4 we have instead for f+that uobeys the equation (u)2=αu2−4ga2u3−2gau4−2gu5/5.(4.58) For this case, the formal solution, y+c=u a dy 1 yα−4ga2y−2gay2−2gy3/5,(4.59) leads to an elliptic function. V. VARIATIONAL APPROACH TO THE FORCED NLSE For the problem of small perturbations to the unforced problem we are interested in variational trial wave functions of the form ψv1(x,t)=A(t)f(βv(t)[x−q(t)])ei{p(t)[x−q(t)]+φ(t)}.(5.1) Here we assume that the collective coordinates (CCs) A(t),φ(t),p(t),q(t) are real functions of time and that fis real. On substituting this trial function in Eq. (2.1) and computing various integrals one finds that the effective Lagrangian is given by L=C0 A2(t) βvp(t)˙ q(t)−˙ φ(t)−[p2(t)−δ]−D1 C0 β2 v +Ck g[A(t)]2(κ+1) βv(κ+1) −4rA(t) cos[kq(t)+φ(t)+θ] ×I[p(t)+k,βv],(5.2) where Ck=∞ −∞ dy[f(y)]2(κ+1),D 1=∞ −∞ [f(y)]2,(5.3) while I[p(t)+k,βv]=∞ 0 dyf (βvy) cos{[p(t)+k]y}.(5.4) Note that for our parametrization of the variational ansatz, ψ∗ψp(t)=1 2i(ψ∗∂xψ−ψ∂xψ∗).(5.5) Thus from our previous discussion about the conservation of momentum in the comoving frame with k=v/2wehavethat p(t) and M(t) are in general not independent but satisfy M(t)[p(t)+k]=const.(5.6) For the case where p(t)+k= 0, one then gets the relation M(t)=const/[p(t)+k].(5.7) A. Variational ansatz Generalizing the κ=1 choice of Ref. [9]forfto arbitrary κ, we choose fv[y]=sech1/κ [y],(5.8) so that our variational ansatz is ψv1(x,t)=A(t) sech1/κ {βv(t)[x−q(t)]}ei{p(t)[x−q(t)]+φ(t)}, (5.9) where A(t) is of the same form as in the unforced solution, namely A(t)=βv(t)2(κ+1) gκ21/(2κ) =β1/κ vα(κ); (5.10) α(κ)=κ+1 gκ21/κ . 046607-7 COOPER, KHARE, QUINTERO, MERTENS, AND SAXENA PHYSICAL REVIEW E 85, 046607 (2012) Note that we have chosen fto be real. In our previous studies of blowup in the NLSE we used a more complicated variational ansatz where there is another variational parameter in the phase multiplying the quadratic term [x−q(t)]2[18,19]. Defining the “mass” of the solitary wave M(t)via M(t)=dxψ∗ψ=C0A(t)2/βv(t)=C0β2/κ−1 vα(κ), (5.11) we then find for the ansatz Eq. (5.10) that for p+k= 0 β2/κ−1 v[p(t)+k]=const.(5.12) For the forcing term contribution to the Lagrangian, Eq. (5.2), we need the integral I[ν,a,βv]=∞ 0 dx cos(ax) sechν(βvx) =[2ν−2 /βv(ν)](ν/2+ia/2βv)(ν/2−ia/2βv), (5.13) where Re(βv)>0, Re(ν)>0, a>0. Special cases of this are obtained for ν=1,2,3. We find I[1,a,βv]= πsechπa 2βv 2βv ,I[2,a,βv]= πacschπa 2βv 2β2 v , I[3,a,βv]=πa2+β2 vsech πa 2b 4β3 v .(5.14) For our variational ansatz Eq. (5.8) C0=√π1 κ 1 2+1 κ;Cκ=√π1+1 κ 3 2+1 κ; (5.15) D1=√π1 κ 2κ23 2+1 κ, so that Cκ=2 2+κC0;D1=1 κ(κ+2)C0.(5.16) Thus we can write everything in terms of C0. B. Arbitrary κ For arbitrary κthe effective Lagrangian is L=M(t)p(t)˙ q(t)−˙ φ(t)−[p2(t)−δ]+β2 v (2 −κ) κ2(2 +κ) −21/κ rβ1/κ−1 v√α(κ)+−cos[kq(t)+φ(t)+θ] 1 κ, (5.17) where ±=1 2±i[k+p(t)] βv(t)+1 κ; (5.18) M(t)=C0α(κ)β2/κ−1. Introducing the notation C=π(k+p) 2β;B=φ+kq +θ, (5.19) an important special solution is obtained in the limit C→0, i.e., p(t)=−k. In this case the Lagrangian becomes L[C=0]=M(t)−k˙ q(t)−˙ φ(t)−(k2−δ)+β2 v (2 −κ) κ2(2 +κ) −21/κ rβ1/κ−1 v√α(κ)21 2κcos(B) 1 κ.(5.20) This leads to the equations ˙ q=−2k, (5.21) ˙ β=−21/κ rα−1/2 κβ−1/κ+1 v1 2+1 κ21 2κsin(B) √π2 κ−121 κ, (5.22) ˙ φ=k2+δ+β2 κ2 −21/κ rα−1/2 κβ−1/κ v(1 −κ)1 2+1 κ21 2κcos(B) (2 −κ)21 κ. (5.23) Assuming the ansatz β=βsand φ(t)=φs−αstin Eqs. (5.22) and (5.23), where βs,αs, and φsare constant, we obtain αs=−2k2,φ s=nπ −kq0−θ, (5.24) where nis an integer and βsis a solution of k2−δ−β2 s κ2 +(−1)n21/κ rα−1/2 κβ−1/κ s(1 −κ)1 2+1 κ21 2κ (2 −κ)21 κ=0. (5.25) These equations become much simpler in the comoving frame, so we will now turn our attention to solving the problem in that frame. VI. VARIATIONAL APPROACH FOR THE ACTION IN THE COMOVING FRAME In the comoving frame we use the following ansatz for the variational trial wave function: vv1(y,t)=A(t)f(βv(t)[y−˜ q(t)])ei{˜ p(t)[y−˜ q(t)]+˜ φ(t)}.(6.1) Comparing with our general ansatz we have the relations y=x+2kt;˜ p=p+k, ˜ q=q+2kt, and ˜ φ=φ+kq +θ. (6.2) On substituting these relations in the effective Lagrangian (5.2), the Lagrangian in the comoving frame takes the form L2=C0 A2(t) βv˜ p(t)˙ ˜ q(t)−˙ ˜ φ(t)−[˜ p2(t)+k2−δ] +β2 v (2 −κ) κ2(2 +κ)+Ck g[A(t)]2(κ+1) βv(κ+1) −4rA(t) cos[ ˜ φ(t)]I[˜ p(t),β].(6.3) 046607-8 FORCED NONLINEAR SCHR ¨ ODINGER EQUATION WITH ... PHYSICAL REVIEW E 85, 046607 (2012) We notice that the canonical momenta to ˜ φand ˜ qare simply given by δL2 δ˙ ˜ q=M(t)˜ p(t); δL2 δ˙ ˜ φ=−M(t).(6.4) A. Variational ansatz function choice Again choosing fv[y] from Eq. (5.8) our variational ansatz for vis vv1(y,t)=A(t)sech1/κ (βv(t)[y−˜ q(t)])ei{˜ p(t)[y−˜ q(t)]+˜ φ(t)}, (6.5) where A(t)isagaingivenbyEq.(5.10). Again the “mass” of the solitary wave M(t) is given by M(t)=dxv∗v= C0A(t)2/βv(t)=C0β2/κ−1α(κ) and from the conservation of Pv=M˜ p[Eq. (2.15)] we find for ˜ p= 0 β2/κ−1˜ p(t)=const.(6.6) We can rewrite Eq. (6.3) as L=M(t)˜ p(t)˙ ˜ q(t)−˙ ˜ φ(t)−[˜ p2(t)+k2−δ]+β2 v (2 −κ) κ2(2 +κ) −21/κ rβ1/κ−1 v√α(κ)+−cos[ ˜ φ(t)] 1 κ,(6.7) where ±=[1 2(±i˜ p(t) βv(t)+1 κ)]. From Lagrange’s equation for ˜ qwe obtain d dt [M(t)˜ p(t)] =0,→β2/κ −1˜ p=β2/κ −1[p(t)+k]=const. (6.8) Letting G(x=˜ p/β)=21/κ r√α(κ)+− 1 κ,(6.9) so that ∂G ∂˜ p=G β,∂G ∂β =−˜ pG β2,(6.10) we get the following differential equations: ˙ ˜ q=2˜ p+1 βM(t)Gβ1/κ−1cos ˜ φ(t),(6.11) M˙ β=Gβ 1/κ κ 2−κsin ˜ φ(t),(6.12) and for ˙ ˜ φwe have 0=(2 −κ) κ˜ p(t)˙ ˜ q(t)−˙ ˜ φ(t)−[˜ p2(t)+k2−δ] +β2 v (2 −κ) κ3+β1/κ M(t)˜ p β2G−1−κ κβ Gcos φ. (6.13) B. ˜ p/β →0 Introducing the notation ˜ C=π˜ p/2β, a special stationary solution is obtained in the limit ˜ C→0, i.e., ˜ p(t)=0. In this case the Lagrangian (6.7) becomes L[˜ C=0] =√π1 κακβv(t)2/κ−1 1 2+1 κ ×δ−˙ ˜ φ(t)−k2+β2 v (2 −κ) κ2(2 +κ) −21/κ r√ακβ1/κ−1 v21 2κcos( ˜ φ) 1 κ.(6.14) This leads to the equations: ˙ ˜ q=0,(6.15) ˙ β=−rAβ(κ−1)/κ 2−κsin( ˜ φ),(6.16) ˙ ˜ φ=δ−k2+β2 κ2−rBβ−1/κ 1−κ 2−κcos( ˜ φ),(6.17) where A=21 κκα−1/2 κ1 2+1 κ21 2κ √π21 κ,A>0, (6.18) B=21 κα−1/2 κ1 2+1 κ21 2κ √π21 κ;B>0. Assuming the ansatz β=βsand ˜ φ(t)=˜ φsin Eqs. (6.16) and (6.17), where βs, and ˜ φsare constant, we obtain ˜ φs=nπ, (6.19) where nis an integer and βsis the solution of k2−δ−β2 s κ2+(−1)nrBβ−1/κ (1 −κ) (2 −κ)=0.(6.20) 1. Linear stability Let us look at the problem of keeping ˜ q=˜ p=0, and looking at small perturbations about the stationary solutions β=βsand ˜ φ=˜ φs=0,π. The analysis depends on whether rcos ˜ φs>0orrcos ˜ φs<0. We discuss the two cases separately. We let ρ=|rcos ˜ φs|>0. Case IA: rcos ˜ φs=ρ>0.The linearized equations for δβ,δφ then become δ˙ β=−Aρβ(κ−1)/κ s 2−κδ˜ φ=−c1δ˜ φ, (6.21) δ˙ ˜ φ=2βs κ2+ρB(1 −κ) κ(2 −κ)β(κ+1)/κ sδβ =c2δβ. (6.22) Combining these two equations by taking one more derivative we find δ¨ β+2δβ =δ¨ ˜ φ+2δ˜ φ=0,(6.23) where 2=c1c2. Whether this will correspond to a stable solution will depend on signs of c1,c2, and hence c1c2. It turns out that the answer depends on the value of κ. In particular, the answer depends on whether κ⩽1or1<κ<2orκ>2. So let us discuss all three cases one by one. Note that the above 046607-9 COOPER, KHARE, QUINTERO, MERTENS, AND SAXENA PHYSICAL REVIEW E 85, 046607 (2012) FIG. 10. Soliton energy computed by subtracting the background energy density from the total energy density. Left panel: r=0(a=0). Right panel: r=−0.075. Parameters are the same as in Fig. 9. where the Lagrangian density reads L=i 2(ψtψ∗−ψ∗ tψ)−|ψx|2+g κ+1(ψψ)κ+1 +δ|ψ|2−re−i(kx+θ)ψ∗−rei(kx+θ)ψ, (7.3) and the dissipation function is given by F=−iα(ψtψ∗−ψ∗ tψ).(7.4) Inserting the ansatz Eq. (5.1) into (7.4) we obtain F=−2αA2f2 v(βv(t)[x−q(t)])[˙ p(q−x)+p˙ q−˙ φ].(7.5) Integrating this expression over space we obtain F=−2αC0 A2 β(p˙ q−˙ φ).(7.6) On the other hand L=dxLis given by Eq. (5.17). Since Fcontains ˙ qand ˙ φ, only the equations for βand pcould be changed by the damping term. For κ=1, we know that the damping only affects the equation for β[8]. For κ=1/2we have the same scenario, i.e., now the equation for βreads ˙ β=−2 3αβ −grC 6β sin(B) sinh(C).(7.7) The factor 2/3 comes from 2/(2/κ −1). VIII. NUMERICAL SIMULATIONS ON THE UNFORCED AND THE FORCED NLSE In this section we would like to accomplish three things. First, we would like to show that our numerical scheme as applied to exact soliton solutions with either r=0orr= 0 leads to known results (r=0) and to results that the exact solutions to the forced problem are metastable and at late times become a solitary wave whose parameters oscillate in time. Second, we want to compare the results of solving the four CC equations for the collective variables β,q,p,φwith a numerical determination of these quantities found by solving the FNLSE. We will find that if ris small compared with the amplitude of the solitary wave, that solution of the CC equations gives a good description of the wave function of the FNLSE. Finally, we would like to demonstrate that the regions of stability for the solitary waves of the FNLSE are well determined by studying the phase portrait as well as the p(v) curve for the approximate solution found by solving the CC equations. Most of the simulations will be restricted to κ=1/2 except for a discussion of blowup of solutions for κ⩾2. A. Numerical methodology Before displaying the results of our simulations, we would like to say a little bit about the numerical method we used and the boundary conditions, since we are constrained to perform the calculation in a box of length 2L.Wehaveallowedthe length to vary from 100 to 400. The number of points used FIG. 11. Time evolution of the exact solution of the FNLSE as in Figs. 9and 10 but using smaller value of the driving terms: k=0.01 and r=−0.005. 046607-16 FORCED NONLINEAR SCHR ¨ ODINGER EQUATION WITH ... PHYSICAL REVIEW E 85, 046607 (2012) FIG. 12. Comparison of the two boundary conditions: periodic (dashed line) vs mixed (solid line) r=−0.075. Left panel: β(t); right panel: soliton profile at t∗=200. The parameters are the same as used in Figs. 9and 10. on the spacial grid was 2L/x. The numerical simulations were performed using a 4th order Runge-Kutta method. N+1 points are used on the spatial grid n=0,1,...N. When studying the exact solutions we have used three different boundary conditions. For r=0 we use “hard wall” boundary conditions, where the wave function vanishes at the boundary: ψ(±L,t)=0.(8.1) For the case of r= 0 and only when studying the exact solutions we have used mixed boundary conditions [see (4.27)]. Otherwise, we have used periodic boundary conditions ψ(−L,t)=ψ(L,t),ψ x(−L,t)=ψx(L,t).(8.2) In one of our simulations we have compared the use of periodic vs mixed boundary conditions and found that the differences in the evolution of both β(t) and ψ(x,t) are hardly visible by the “eyeball method” (see Fig. 12) The other parameters related with the discretization of the system are increasing values of L, namely L=50, L=62.8, L=100,or L=200. We have chosen as our grids in xand t,x =0.05 and t =0.0001 [such that t < (1/2)(x)2]. We next need to determine the collective variables qand βused in the CC equations. We determine q(t) from our numerical simulation by equating it to the value of xfor which the density of the norm |ψ|2is maximum. To determine β(t) for finite energy solitary waves, we assume that the variational parametrization of ψ(x,t), namely Eq. (5.9) with Agiven by Eq. (5.10) is an adequate description of the wave function. If we do this and determine ψ[x=q(t)] then we have that |ψ|2[x=q(t)] =A2=β2(κ+1) gκ21/κ .(8.3) Then we determine βas follows: β={gκ2|ψ|2κ[x=q(t)]}/(κ+1).(8.4) In the particular case when we are studying the time evolution of the exact solution for κ=1/2, i.e., Eq. (4.16), with conditions of Eq. (4.20), we need to subtract off the constant term to determine β(t). From Eq. (4.16) we can obtain β(t) from |ψ|2|x=q(t)=a+6β(t)2 g.(8.5) One can also compute the momentum P(t) given by Eq. (2.12). As an initial condition in most of our simulations (when we are not discussing the exact solution), we will use an approximate solution given by the variational ansatz Eq. (5.9) with initial conditions for β0,p0,q0=0, and φ0. In comparing with the CC equations, we will solve the four ordinary differential equations for β,p,q,φwhich for arbitrary κ are given in Eqs. (6.8)–(6.13). FIG. 13. Here we compare the time evolution of a solitary wave traveling to the left initially with constant velocity −2kstarting from an approximate stable stationary solution for three different values of positive r. The times are t∗=33.3, 66.6, 100, 133.3, 166.6, 200. We use mixed boundary conditions. Left panel: r=0.01; middle panel: r=0.05; right panel: r=0.075. Parameters: κ=1/2, δ=−1, g=2, k=0.1. Initial condition (4.16) with (4.19). 046607-17 COOPER, KHARE, QUINTERO, MERTENS, AND SAXENA PHYSICAL REVIEW E 85, 046607 (2012) FIG. 14. Time evolution of the inverse width parameter β(t). Here we plot β(t) for the same initial conditions as in Fig. 13.Againwehave left panel r=0.01, middle panel r=0.05, and right panel r=0.075. B. Numerical simulations of PDE for r=0, arbitrary κ The initial conditions are those of the exact one-soliton solution of the unforced NLSE given by Eq. (3.1).Forr=0 and δ=0, the stability of the NLSE has been well studied as a function of κ.Forκ<2 the solutions are known to be stable, and for κ>2 the solutions are unstable. For κ=2 there is a critical mass above which the solutions are unstable. The nature of the solution when it is unstable has been studied in variational approximations [19] as well as using various numerical algorithms [23–26]. Here we want to show that our code reproduces these well known facts. We are also interested in the effect of δin our simulations and we find that the critical value of κis independent of δ. This agrees with our arguments earlier on the effect of scale transformations on the stability of the exact solutions for r=0. That is, stability does not depend on δ. C. Simulations at r=0, different values of κ In this section we study the numerical stability of the exact solutions for r=0 for different values of κ(κ∈[0.25,3]). FIG. 15. (Color online) Test of stable stationary solution of the CC equations (horizontal line) by a simulation of the NLSE with periodic boundary conditions (black and red-upper solid lines for system sizes 125.6 and 251.2, respectively, the latter line is shifted upwardby0.05). The influence of linear excitations, which are radiated at early times, on the soliton oscillations is discussed in the text. The parameters used are κ=1/2, r=0.05, k=−0.1, θ=0, g=2, δ=−1, α=0, with initial conditions p0=−k,q0=0, φ0=0, and β0=βs=0.517 73. For the CC equations we choose p0=−k+10−5to avoid numerical singularities. The approximate solitary wave is described by Eq. (6.69), and is shown in Fig. 5. First we set r=0, δ=0, g=2 in NLSE and we start from the exact soliton solution (3.1) with β=1, p=0.01 and φ0= 0. In this simulation the boundary conditions chosen were (±L,t)=0. We notice for κ<2 as shown in Fig. 6the solitary wave moves to the right and maintains its shape. Indeed, in the simulations the solitons are not stable for κ⩾2 (see Fig. 7), the amplitude of the unstable soliton grows and the soliton becomes narrow. Also the soliton moves only slowly to the right while blowing up at a finite time. In Fig. 8, upper panel, we show the real and imaginary parts of the field for κ=1.5,2,3 for the final time of the simulations. The result of adding a term proportional to δdoes not alter the behavior of the solitary waves. Choosing δ=−1 we obtain similar results as for δ=0, i.e., the soliton is unstable for κ⩾2. In Fig. 8, right lower panel, we show the evolution of βfor short times and κ=1.5; 2; 3. This shows that once we reach the metastable solution for κ=2, we see β(t)→∞in a finite amount of time. This agrees with our analysis based on Derrick’s theorem [see Eq. (3.12)], and with the discussion of the critical mass needed for blowup for κ=2inRef.[19]. For δ=−1wehave also investigated the constant phase solutions that fulfill the condition (3.5). We have studied numerically the case r=0, δ=−1, g=2 with initial conditions p=0.01, φ0=0, and β=−κ2(δ+p2) and varied κ∈[0.25,3]. We found that again the soliton becomes unstable for κ⩾2 in accord with our result that instability as a function of κdoes not depend on δ. Our numerical experiments for the case r=0 show that our codes reproduce well known results for the stability of the solutions. In a future paper we will compare the numerical solutions in the unstable regime with the predictions of the CC method. Here we will use the exact form of the solution rather than the post-Gaussian trial functions we used earlier [19] as well as include an additional variational parameter in the phase that is canonically conjugate to the width parameter β. To estimate the finite size effect on the definition of βwe notice in the left panel of Fig. 9that the values of β(t)inthe simulation of the NLSE deviate from the exact constant value of 0.5 by approximately 0.25%. By increasing Lone could reduce this error. D. On the exact solutions for r= 0andκ=1/2 Earlier we showed that for r= 0 and κ=1/2 there are exact solutions of the form [see Eq. (4.24)] f(y)=−sgn(r)(a+bsech2βy),(8.6) 046607-18 FORCED NONLINEAR SCHR ¨ ODINGER EQUATION WITH ... PHYSICAL REVIEW E 85, 046607 (2012) FIG. 16. Comparison of simulations of the collective coordinate (CC) approximation and the numerical solution of the PDE where we have considered the case κ=1/2 and included damping when appropriate. Left upper panel: Soliton moving to the left for t∗=250,500. Right upper panel: Real (solid line) and imaginary (dashed line) parts of ψ(x,t)fort∗=500 (both lines are shifted down by 0.5). Real (solid line) and imaginary (dashed line) parts of the exact solution for the undamped NLSE for t=500 (both lines are shifted upward by 0.5). Left and right middle panels: comparison of the time evolution of βand qcomputed from the simulations of the PDE (solid line) and numerical solutions of the CC equations (dashed line). Lower panels (different scales are shown): for t∗=500 soliton profile from simulations (solid line) and from the exact solution (4.16) and (4.19) (dashed line). Notice that the exact solution is obtained for α=0. Parameters: g=2, k=0.1, r=0.05, δ=−1, θ=0, α=0.05, with initial conditions β0=1, p0=0, q0=0, and φ0=−1.69 +π/2≈−0.119. Notice that for this set of parameters the condition |r|<k 4/(4g)[seebelowEq.(4.19)] is satisfied. After a transient time, since the numerical solution approaches the stationary solution, the p(v) curve is just a point (−0.1,−0.2). The phase portrait is also represented by a point (0.044 911 697 24,−0.028 190 023 64). with a,β, and bgiven by Eq. (4.20). In the numerical simulation shown in Figs. 9and 10 the parameters used in the simulation are κ=1/2,δ =−1,g =2,k =0.1,r =−0.075. For that case, f(y)=0.486 113 sech2(0.402 539y)+0.090 462 2.(8.7) The unstable stationary variational solution corresponding to this exact solution is f2+(y)=0.745 535 sech2(0.498 509y).(8.8) These solutions are shown in Fig. 2. As we can see from the numerical solution of the FNLSE shown in Fig. 9, the exact solitary wave oscillates in amplitude and width which shows that the original exact solution was unstable. In spite of this the solitary wave neither dissipates nor blows up. The width parameter oscillates from its initial value of around 0.4toan upper value of around 0.8 with a period Tof around T=20. If one reduces the driving terms to k=0.01 and r=−0.005, the solitary wave again oscillates in amplitude and width but the amplitude of the oscillation is only 10% of the total amplitude 046607-19 COOPER, KHARE, QUINTERO, MERTENS, AND SAXENA PHYSICAL REVIEW E 85, 046607 (2012) FIG. 17. Results from simulating the time evolution of the solitary wave using both the PDE as well as the CC equations for κ=1/2 with no damping. We use the same parameters as in Fig. 16,butα=0. Upper panels: soliton profiles (different scales) for t=250 (solid line) and t=500 (dotted line). Middle panels: time evolution of βcomputed from the simulations of PDE (left) and computed from the numerical solutions of CC equations (right). Lower left panel: elliptic orbit in the phase portrait with positive sense of rotation which predicts stability; lower right panel: positive slope of the p(v) curve predicts stability. and the oscillation frequency is reduced from the previous case as is shown in Fig. 11. In Fig. 12 we compare the behavior of β(t), and the late form of the wave functions using two types of boundary conditions (BCs): periodic BCs [ψ(−L,t)=ψ(L,t), ψx(−L,t)= ψx(L,t)] shown as dashed lines, and “mixed” BCs which are shown as solid lines. In this simulation the two boundary conditions lead to almost identical results. Most of the simulations (except the ones starting from the exact solution) were performed using periodic boundary conditions. Next we look at the case where r>0 where for amplitudes A(t)>0 there are no exact solutions. However, the stationary solutions of the CC equations that are stable to linear perturbations do lead to solitary waves that have widths whose oscillations have only very small amplitude. To show that we start with the stationary solutions with different positive rbut other parameters (κ,δ,g,k) the same as in Fig. 9. We show this in Fig. 13 where we consider solitary waves moving to the left initially with constant velocity v=−2k. We see that as we increase rthe amplitude of the solitary wave increases and in Fig. 14 we see that the average value of βalso increases with r.WenoticeinFig.14 that the oscillations in βget more irregular as we increase the value of r. E. Comparison of numerical simulations of the FNLSE with the solution of the equations for the collective coordinates In this section we would like to show that the solution to the CC equations gives quantitatively good results for the collective variables p(t),β(t),q(t), and φ(t) when these quantities are calculated directly from the solution of the FNLSE. We also want to show that the criterion for stability of the CC equations, namely a study of the phase portrait or the p(v) curve, gives an accurate measure of what initial 046607-20 FORCED NONLINEAR SCHR ¨ ODINGER EQUATION WITH ... PHYSICAL REVIEW E 85, 046607 (2012) FIG. 18. Unstable stationary solution for r=0.05, κ=1/2, g=2, k=−0.1, θ=0, δ=−1, and α=0. The stationary solution has p0=−k(simulations) and p0=−k+10−5(numerical solutions of the CCs), φ0=π,andβ0=βs=[(k2+k4−8rg/3)/8]1/2[see Eq. (6.70)]. Here we find for β(t) that both the simulation (solid line) and the solution to the CC equations (dotted line) show switching to a new solution. From the p(v) curve the negative slope of the curve predicts this instability. In addition, the orbit of the phase portrait is a separatrix, which also predicts the instability seen in β. values of collective variables lead to stable vs unstable solitary wave solutions. First we discuss the linearly stable stationary solution Eq. (6.69) that is shown in Fig. 5. In the CC equation βremains at the fixed value β(0) =βs=0.517 73. Using this as an initial condition in the FNLSE, solved using periodic boundary conditions, one finds that the solitary wave has a mean value of βonly 2% different from the solution, Fig. 15. For 0 <t⩽30, phonons (short for linear excitations) are radiated by the soliton which travel faster than the soliton. Coming from a boundary these phonons reappear and collide with the soliton producing the increased oscillations seen in Fig. 15. For the system of length 2L=125.6, these increased oscillations begin at tc≈125. When the system size is doubled, both the soliton and the phonons have to cover doubled distances before the collisions begin at tc≈250. If we add a little bit of damping and start from a value of β0which dissipates to a stationary point then all the collective variables are well approximated by the solution of the CC equation as shown in the middle panel of Fig. 16.Herewehave chosen g=2, k=0.1, r=0.05, δ=−1, θ=0, α=0.05. For these values, the stationary solution is approximately given by Eq. (6.68) which has a value of βs=0.505 12. In the middle panel of Fig. 16 we use the initial condition β0=1.We find that both the CC equations and the solution of the FNLSE relax to the approximate stationary solution of Eq. (6.68). In Fig. 17 we turn off the damping and see how both the simulation and CC equations evolve. Solving the CC equations with the initial condition β0=1, one finds from both the orbit of this initial condition in the phase portrait (see left bottom panel) and the p(v) curve that the solitary wave should be stable. In the middle panel we see that both the exact solution and the solution to the CC equation for βoscillate with 1% oscillations around either the initial value β0=1(forthe numerical solution) and a slightly shifted value (0.994) for the CC equation. The actual solution has another oscillation frequency for the height at maximum. In Fig. 18 we study the time evolution of a stationary solution that is known to be unstable. We use the solution f2+(y)giveninEq.(6.70). Both the CC equations and the simulation of the FNLSE show that this solitary wave develops into a solution that is represented by a separatrix in the phase portrait. The negative slope of the p(v) curve also predicts instability. The soliton stability depends strongly on β0.Forthe parameters and initial conditions of Fig. 19 both stability criteria predict the following pattern: instability for β0⩽ 0.484 511, stability for 0.484 511 <β 0⩽0.5458, instability for 0.5458 <β 0⩽0.65, and stability for β0⩾0.66. In Fig. 19 we present two examples: a stable soliton (β0=0.53) and an unstable one (β0=0.65), both confirmed by our simulations. We note, however, that the above boundaries between stable and unstable regions do not always fully agree with the simulations: the errors vary between 1.7% and 12%. In Fig. 20 we show a particular case where the time evolution of an initial condition which was close to an unstable 046607-21 COOPER, KHARE, QUINTERO, MERTENS, AND SAXENA PHYSICAL REVIEW E 85, 046607 (2012) FIG. 19. (Color online) Simulations for the NLSE with r=0.05, k=−0.1, periodic BC, with initial conditions β0=0.53 (black curves) and 0.65 (red, or lower solid grey curve). The other parameters and initial conditions are κ=1/2, g=2, δ=−1, θ=0, α=0, φ0=0, p0=−k+0.001,and q0=0. Upper left panel: β(t) compared with numerical solutions of the CC equations (dotted lines). Upper right panel: orbits of the phase portrait (red horseshoe for β0=0.65 and black ellipse for β0=0.53). Lower left panel: p(v)curveforβ0=0.53. Lower right panel: p(v)curveforβ0=0.65. stationary solution of the CC equation exhibits intermittency in both the inverse width (and thus the amplitude Aalso) as well as the energy. We have also observed intermittency in the solutions of the CC equations for κ=1. IX. CONCLUSIONS In this paper we have studied analytically, numerically, and in a variational approximation the FNLSE with arbitrary nonlinearity parameter κ. We studied in detail a variational approximation based on the solutions to the unforced problem and have studied the CC equations coming from the Euler- Lagrange equations. We determined the stationary solutions of the CC equations and studied their linear stability. We found that for small forcing parameter, the CC equations give quantitative agreement with directly solving the FNLSE. Also the domains of stability of the initial conditions for the FNLSE were quite remarkably close to those found by studying the orbits in the phase portrait for the CC equations as well as the stability curves p(v). We also found that the linearly stable stationary solutions of the CC equations were quite close to the stationary solutions of the FNLSE, and when these stationary solutions were used as initial conditions for the FNLSE the oscillations of the width parameter were of small amplitude. These simulations at κ=1/2 reinforce our belief that our two dynamical criteria for understanding the stability of the solitary wave, namely the stability curve p(v) as well as studying orbits in the phase portrait are accurate indications of the stability of the solitary waves as obtained by numerical simulations of the FNLSE. Our results are likely to shed light on the behavior of a number of physical systems varying from optical fibers [1], Bragg gratings [2], BECs [3,4], nonlinear optics, and photonic crystals [7]. We intend to extend our study of the FNLSE equation to the regime κ⩾2 to see how forcing affects the blowup of solitary wave solutions in that regime. This will be done in the variational approximation by adding a variational parameter canonically conjugate to the width parameter β similar to the approach of Refs. [18,19]. In this paper we have considered only the simplest forcing term which is purely harmonic. More elaborate forcing terms such as nonparametric spatial-temporal driving forces are possible, which have been discussed for the special case κ=1inRef.[8]. We have conjectured the dynamic criteria for the standard NLSE and used it for the generalized NLSE with arbitrary nonlinearity (κ)for determining the stability of solitary waves: the p(v) curve and the properties of the orbits in the phase portrait of soliton. We are currently exploring whether this approach can be applied to other forced nonlinear problems. We have considered FNLSE for the one-dimensional case. It is conceivable that our analysis may hold for quasi-one-dimensional solitary wave solutions of the corresponding three-dimensional equation. Also, we 046607-22 FORCED NONLINEAR SCHR ¨ ODINGER EQUATION WITH ... PHYSICAL REVIEW E 85, 046607 (2012) FIG. 20. Numerical simulation of the NLSE for r=0.05 using periodic boundary conditions and κ=1/2. β(t)andE(t) show intermittency. Upper panels: Soliton moving to the right for t∗=1000 and t∗=2000 and the evolution of the energy of the system. Lower panels: time evolution of βand qcomputed from the simulations of PDE. Other parameters of the simulations: g=2, k=−0.01, δ=−1, θ=0, α=0, with initial conditions p0=ps=−k,q0=0, φ0=π,andβ0=βs=0.4983. This corresponds to the unstable stationary solution of Eq. (6.72), namely f2+(y)=0.745 042 sech2(0.498 345y). did not consider the uniqueness of the soliton solutions here since it is highly nontrivial to assess uniqueness for nonlinear equations. All these points merit further study. ACKNOWLEDGMENTS This work was supported in part by the US Department of Energy. F.G.M. acknowledges the hospitality of the Mathematical Institute of the University of Seville (IMUS) and of the Theoretical Division and Center for Nonlinear Studies at Los Alamos National Laboratory and financial support by the Plan Propio of the University of Seville and by Junta de Andaluc´ ıa. N.R.Q. acknowledges financial support by the MICINN through FIS2011-24540, and by Junta de Andaluc´ ıa under Projects No. FQM207, No. FQM-00481, No. P06-FQM-01735, and No. P09-FQM-4643. [1] C. De Angelis, IEEE J. Quantum Electron. 30, 818 (1994). [2] J. Atai and B. A. Malomed, Phys. Lett. A 284, 246 (2001). [3] F. Kh. Abdullaev, A. Gammal, L. Tomio, and T. Frederico, Phys. Rev. A 63, 043604 (2001). [4] W. Zhang, E. M. Wright, H. Pu, and P. Meystre, Phys. Rev. A 68, 023605 (2003). [5] A. Minguzzi, P. Vignolo, M. L. Chiofalo, and M. P. Tosi, Phys. Rev. A 64, 033605 (2001); B. Damski, J. Phys. B 37, L85 (2004). [6] F. Kh. Abdullaev and M. Salerno, Phys.Rev.A72, 033617 (2005). [7] I. Maruno, Y. Ohta, and N. Joshi, Phys. Lett. A 311, 214 (2003). [8] F. G. Mertens, N. R. Quintero, and A. R. Bishop, Phys. Rev. E 81, 016608 (2010). [9] F. G. Mertens, N. R. Quintero, I. V. Barashenkov, and A. R. Bishop, Phys.Rev.E84, 026614 (2011). [10] N. R. Quintero, F. G. Mertens, and A. R. Bishop, Phys. Rev. E 82, 016606 (2010). [11] U. Peschel, O. Egorov, and F. Lederer, Opt. Lett. 15, 1909 (2004). [12] D. J. Kaup and A. C. Newell, Proc. R. Soc. London A 361, 413 (1978); Phys. Rev. B 18, 5162 (1978). [13] P. S. Lomdahl and M. R. Samuelsen, Phys.Rev.A34, 664 (1986). [14] B. A. Malomed, Phys. Rev. E 51, R864 (1995); G. Cohen, ibid. 61, 874 (2000). [15] K. Nozaki and N. Bekki, Physica D 21, 381 (1986). [16] G. H. Derrick, J. Math. Phys. 5, 1252 (1964). [17] R. Jackiw and A. Kerman, Phys. Lett. A 71, 158 (1979); F. Cooper, S. Y. Pi, and P. N. Stancioff, Phys.Rev.D34, 3831 (1986). [18] F. Cooper, C. Lucheroni, and H. K. Shepard, Phys. Lett. A 170, 184 (1992). [19] F. Cooper, C. Lucheroni, H. K. Shepard, and P. Sodano, Physica D68, 344 (1993). 046607-23 COOPER, KHARE, QUINTERO, MERTENS, AND SAXENA PHYSICAL REVIEW E 85, 046607 (2012) [20] F. Cooper, A. Khare, B. Mihaila, and A. Saxena, Phys. Rev. E 82, 036604 (2010). [21] I. V. Barashenkov, T. Zhanlav, and M. M. Bogdan, in Nonlinear World. IV International Workshop on Nonlinear and Turbulent Processes in Physics,editedbyV.G.Bar’yakhtar, M. Chernousenko, N. S. Erokhin, A. G. Sitenko, and V. E. Zakharov (World Scientific, Singapore, 1990), pp. 3–9. [22] I. V. Barashenkov and Yu. S. Smirnov, Phys. Rev. E 54, 5707 (1996). [23] H. A. Rose and M. I. Weinstein, Physica D 30, 207 (1988). [24] T. C. Scott, R. A. Moore, G. J. Fee, M. B. Monagan, and E. R. Vrscay, J. Comput. Phys. 87, 366 (1990). [25] W. Kutzelnigg, Z. Phys. D 15, 27 (1990); in Relativistic Electronic Structure Theory, Part I, edited by P. Schwerdtfeger (Elsevier, Amsterdam, 2002). [26] C. I. Siettos, I. G. Kevrekidis, and P. G. Kevrekidis, Nonlinearity 16, 497 (2003). 046607-24