A Useful integrals and definitions

# Stability of new exact solutions of the nonlinear Schrodinger equation in a Poschl-Teller external potential

## Abstract

We discuss the stability properties of the solutions of the general nonlinear Schrödinger equation (NLSE) in 1+1 dimensions in an external potential derivable from a parity-time () symmetric superpotential that we considered earlier Kevrekidis et al. (2015). In particular we consider the nonlinear partial differential equation for arbitrary nonlinearity parameter , where and is the well known Pöschl-Teller potential which we allow to be repulsive as well as attractive. Using energy landscape methods, linear stability analysis as well as a time dependent variational approximation, we derive consistent analytic results for the domains of instability of these new exact solutions as a function of the strength of the external potential and . For the repulsive potential (and ) we show that there is a translational instability which can be understood in terms of the energy landscape as a function of a stretching parameter and a translation parameter being a saddle near the exact solution. In this case, numerical simulations show that if we start with the exact solution, the initial wave function breaks into two pieces traveling in opposite directions. If we explore the slightly perturbed solution situations, a 1% change in initial conditions can change significantly the details of how the wave function breaks into two separate pieces. For the attractive potential (and ), changing the initial conditions by 1 % modifies the domain of stability only slightly. For the case of the attractive potential and negative perturbed solutions merely oscillate with the oscillation frequencies predicted by the variational approximation.

1

## I Introduction

The study of open systems with balanced loss and gain, typically defined by Parity-Time () symmetry, has elicited significant attention from physics, nonlinear science and mathematics communities during the past decade. This is in part due to their emerging applications in many physical contexts and in part due to their intriguing mathematical structure. The initial investigation of such systems Bender (2007); Geyer et al. (2006); Fring et al. (2008); Bender et al. (2012) arose in the context of whether non-Hermitician quantum systems could lead to entirely real eigenvalues. Keeping in perspective the formal similarity of the Schrödinger equation with Maxwell’s equations in the paraxial approximation, many experimentalists realized that such -invariant systems can indeed be fabricated using optical means Makris et al. (2011); Ruschhaupt et al. (2005); Makris et al. (2008); Klaiman et al. (2008); Longhi (2009a, b, 2010); Ruter et al. (2010); Guo et al. (2009); Regensburger et al. (2012). Motivated by this success, in the ensuing years, -invariant phenomena was also observed in electronic circuits Schindler et al. (2011, 2012), mechanical constructs Bender et al. (2013), whispering-gallery microcavities Peng et al. (2014), among many other physical systems.

In a parallel development, the concept of supersymmetry (SUSY) prevalent in high-energy physics was also experimentally studied in optics Miri et al. (2013); Heinrich et al. (2014). The underlying notion is that for a given potential we can obtain a SUSY partner potential such that both potentials possess identical spectrum (with possibly one eigenvalue different) Gendenshteĭn and Krive (1985); Cooper et al. (1995). A simultaneous presence of -symmetry and SUSY can lead to unexpected phenomena and is likely to be very useful in achieving transparent and one-way reflectionless complex optical potentials Bagchi and Roychoudhury (2000); Bagchi et al. (2001); Bagchi and Quesne (2000); Ahmed (2001); Midya (2014). Previously Kevrekidis et al. (2015) we studied the interplay between nonlinearity, -symmetry and supersymmetry as well as the rich consequences of this interplay. There we obtained exact solutions of the general nonlinear Schrödinger equation (NLSE) in 1+1 dimensions in the presence of a -symmetric complex potential Gendenshteĭn and Krive (1985); Cooper et al. (2001). In a recent paper Cooper et al. (2017) we studied the stability properties of the solutions of NLSE in the real partner potential of the problem studied in Kevrekidis et al. (2015) which was a Pöschl-Teller potential Poschl and Teller (1933); Landau and Lifshitz (1989).

Here our objective is to discuss the stability properties of two related exact solutions which exist when we change the sign of the nonlinear coupling to being negative keeping the potential attractive, or keep the sign of the nonlinear term unchanged but consider a repulsive potential. In the latter case we will find that the solutions are translationally unstable, whereas in the former case the solutions are stable to small perturbations

### i.1 Different solutions to the NLSE in an external Pöschl-Teller potential

By allowing the nonlinearity coupling and the sign of the potential we have found different classes of exact soltuions when the NLSE is in the presence of a Pöschl-Teller potential centered at . Schrödinger’s equation for these cases is given by:

 {i∂t+∂2x+g\absψ(x,t)2κ−V(x)}ψ(x,t)=0, (1)

where

 V(x)=−λ~b2\sech2(x)\qc~b2=b2−1/4, (2)

with and . Since can be scaled out of the equation by letting , we can restrict ourselves to in what follows. The signs here are chosen such that the nonlinear term is attractive for and repulsive for and the external Pöschl-Teller potential is attractive for and repulsive for . This potential is a special case of potentials obtainable from the complex -symmetric SUSY superpotential

 W(x)=\qty(m−1/2)tanhx−ib\sechx, (3)

with , which gives rise to -symmetric partner potentials . Our real corresponds to . There are several cases of Eq. (1) which have exact solutions. These are

1. Attractive nonlinear term and attractive potential: , . In this case, the exact solution is given by

 ψ0(x,t) =A0(~b,γ)\sechγ(x)eiγ2t (4a) A2/γ0(~b,γ) =γ(γ+1)−~b2, (4b)

where . In this case,

 ~b2γ≡γ(γ+1)≥~b2≥0. (5)

We studied this case in a previous paper Cooper et al. (2017), where we found that all solitary waves for and are stable, as for the case of solitary waves in the NLSE (). However, we also found a new region above where these solutions are stable.

2. Attractive nonlinear term and repulsive potential: , . For this case, the exact solution is given by

 ψ0(x,t) =A0(~b,γ)\sechγ(x)eiγ2t (6a) A2/γ0(~b,γ) =γ(γ+1)+~b2. (6b)

In this case, we only require . This solution goes over to a particular moving solitary wave solution of the NLSE when . Since the solutions of the NLSE are stable to deformations of the width for all we expect (and we will find) that in that regime there will be a critical value of above which the solution will be unstable to width deformations. We expect and we find that again the solutions are always unstable for . What we will also find is that for all values of these solutions are unstable to a slight translation, even if induced by numerical noise.

3. Repulsive nonlinear term and attractive potential: , . For this case, the exact solution is given by

 ψ0(x,t) =A0(~b,γ)\sechγ(x)eiγ2t (7a) A2/γ0(~b,γ) =~b2−γ(γ+1). (7b)

In this case, we require . For this choice of there are no solitary wave solutions in the absence of the potential. We will find that these solutions are linearly stable.

In all these cases, we find that the quantity

 gA2/γ0(~b,γ)=γ(γ+1)−λ~b2, (8)

is independent of and depends only on the sign of . The normalization, or “mass” of these exact wave functions is given by

 M0(~b,γ)=∫\ddx\absψ0(x,t)2=A20(~b,γ)c1[γ], (9)

where

 c1[γ]=∫\ddz\sech2γ(z)=√πΓ[γ]Γ[γ+1/2]. (10)

This paper is structured as follows: in Section II we discuss Hamilton’s principle of least action and the time-dependent variational approximation. In Section III, we use Derrick’s theorem to study the stability of these solutions to width instabilities. In Section IV we discuss the energy landscape when we include translations of the origin of the solution. In Section V we perform a linear stability analysis. In Section VI, we introduce a four-parameter trial wave function to study the dynamics of the model, and in Section VII we provide results of the direct numerical solutions of the nonlinear Schrödinger equation in the Pöschl-Teller external potential. Our main conclusions are summarized in Section VIII.

## Ii Time-dependent variational principle

The time-dependent version of the variational approximation can be traced to an obscure appendix in the 1930 Russian edition of the “Principles of Wave Mechanics,” by Dirac.2 In this version of the variational approximation, the wave function is taken to be a function of a number of time-dependent parameters. Variation of the action, as defined by Dirac, leads to a classical set of Hamiltonian equations of motion for the parameters. These classical equations are then solved as a function of time to provide an approximation to the evolution of the wave function.

The action which leads to Eq. (47) is given by

 Γ[ψ,ψ∗]=∫\ddtL[ψ,ψ∗] (11)

where

 L[ψ,ψ∗] =i2∫\ddx\qty[ψ∗(∂tψ)−(∂tψ∗)ψ]−H[ψ,ψ∗], (12a) H[ψ,ψ∗] =∫\ddx[\abs∂xψ2−g\absψ2κ+2κ+1+V(x)\absψ2]. (12b)

The NLSE and its complex conjugate follow from minimizing the action via,

 δΓδψ∗=δΓδψ=0. (13)

### ii.1 Symplectic formulation

In this section it will be useful to introduce a symplectic formulation of Lagrange’s equations for the variational parameters. We consider a variational wave function of the form:

 ~ψ[x,Q(t)]\qcQ(t)={Q1(t),Q2(t),…,Q2n(t)}.

Introducing the notation , the Lagrangian (12a) is given by

 L[Q,˙Q]=πi(Q)˙Qi−H[Q], (14)

where

 πi(Q)=i2∫\ddx{~ψ∗[∂i~ψ]−[∂i~ψ∗]~ψ}, (15)

and is given by

 H(Q)=∫\ddx[|∂x~ψ|2−g|~ψ|2κ+2κ+1+V(x)|~ψ|2]. (16)

The Euler-Lagrange equations now become

 \dvt(\pdvL˙Qi)−\pdvLQi=0. (17)

From (14) this gives

 fij(Q)˙Qj=vi(Q), (18)

where we have set , and where

 fij(Q)=∂iπj(Q)−∂jπi(Q) (19)

is an antisymmetric symplectic matrix. If , we can define an inverse as the contra-variant matrix with upper indices,

 fij(Q)fjk(Q)=δik, (20)

in which case the equations of motion (18) can be put in the form:

 ˙Qi=fij(Q)vj(Q). (21)

Conservation of energy is expressed as

 \dvH(Q)t=˙Qivi(Q)=fij(Q)vj(Q)vi(Q)=0, (22)

since is an antisymmetric tensor. Poisson brackets are defined using this antisymmetric tensor. If and are functions of , Poisson brackets are defined by

 {A(Q),B(Q)}=(∂iA(Q))fij(Q)(∂jB(Q)). (23)

In particular,

 {Qi,Qj}=fij(Q). (24)

This definition satisfies Jacobi’s identity. That is, what we have shown here is that the quantities are classical symplectic variables.

## Iii Derrick’s theorem

Derrick’s theorem Derrick (1964) states that for a Hamiltonian dynamical system, an exact solution of the equation of motion is unstable if under a scale transformation, with fixed normalization, the energy of the system is lowered. The stretched wave function for Derrick’s theorem is given by

 ψβ(x,t)=A(~b,β,γ)\sechγ(βx), (25)

with the normalization fixed by the requirement,

 M[~b,β,γ] =∫\ddx\absψβ(x,t)2=A2(~b,β,γ)c1[γ]/β (26) =M0[~b,γ]=A20(~b,γ)c1[γ].

So . Evaluation of the Hamiltonian (12b) with Derrick’s wave function gives:

 H(β,γ)=H1(β,γ)+H2(β,γ)+H3(β,γ), (27)

where

 H1(β,γ) =∫\ddx\abs∂xψβ2 (28a) =A20β2γ2∫\ddz\sech2γ+2(z)sinh2(z) =A20β2γ2c1[γ+1], H2(β,γ) =−gκ+1∫\ddx\absψβ2κ+2 (28b) =−gA2/γ0A20γβ1/γγ+1∫\ddz\sech2γ+2(z), =−A20γβ1/γγ+1[γ(γ+1)−λ~b2]c1[γ+1], H3(β,γ) =∫\ddxV(x)\absψβ2 (28c) =−λ~b2A20β∫\ddx\sech2(x)\sech2γ(βx) =−λ~b2A20g1[β,γ],

where

 g1[β,γ]=∫\ddz\sech2γ(z)\sech2(z/β). (29)

So then is given by

 h(~b,β,γ)≡H(β,γ)A20(~b,γ)c1[γ+1] (30) =12β2γ−γβ1/γγ+1[γ(γ+1)−λ~b2]−λ~b2g1[β,γ]c1[γ+1],

and is independent of . is stationary when . We have

 \pdvh(~b,β,γ)β=γ[β−β1/γ−1] (31) +λ~b2\qty[β1/γ−1γ+1−1c1[γ+1]\pdvg1[β,γ]β].

From Eq. (98) in Appendix A, we find

 \pdvg1[β,γ]β∣∣β=1=c1[γ+1]γ+1, (32)

so that at ,

 \pdvh(~b,β,γ)β∣∣β=1=0, (33)

for all values of and . The sign of the second derivative of with respect to at determines whether the solution is unstable to small changes in the width. If evaluated at is negative the solution is unstable. We find

 \pdv[2]h(~b,β,γ)β=γ+(γ−1)β1/γ−2 (34) −λ~b2\qty[γ−1γ(γ+1)β1/γ−2+1c1[γ+1]\pdv[2]g1[β,γ]β].

From Appendix A, we have

 \pdv[2]g1[β,γ]β∣∣β=1=4c2[γ+1]−6c2[γ+2]−2c1[γ+1]γ+1,

where

 c2[γ] =∫\ddzz2\sech2γ(z) (35) =22γ−14F3[γ,γ,γ,2γ;1+γ,1+γ,1+γ;−1]/γ3.

Inserting this into (34) and evaluating it at gives:

 \pdv[2]h(~b,β,γ)β∣∣β=1=2γ−1 (36) +λ~b2\qty[1γ−2γ+1γ2c2[γ+1]−3c2[γ+2]c1[γ]],

which is independent of .

The critical value of , when (36) vanishes, is

 ~b2crit=−λγ(2γ−1)1−(2γ+1)2c2[γ+1]−3c2[γ+2]c1[γ], (37)

which is independent of . One can easily check that

 1−(2γ+1)2c2[γ+1]−3c2[γ+2]c1[γ]>0, (38)

for all . In Appendix A, we give an alternative form for Eq. (37), which is in agreement with Ref. [Cooper et al., 2017]. In Fig. b, we plot for the three cases. Referring to the figure, according to Derrick’s theorem,

1. For case I with , , and , we see from Eq. (36) that for all , so Derrick’s theorem predicts that solutions are width stable for . For , there is another region for where width stable solutions are also possible.

2. For case II with , , and , width stable solutions are possible for , but we will find that this region is unstable to translation perturbations.

3. For case III with , , and width stable solutions are possible for all .

Derrick’s theorem is a version of the time-independent variational method and only provides information about the stability of the system under a change in , the width of the wave function. Thus Derrick’s theorem only gives a sufficient condition for an instability to occur. To see if there are translational instabilities as well as width instabilities we will consider what happens to the energy when we make a small translation of the position of the solution.

## Iv Translational Stability Landscape

Using Derrick’s theorem we explored whether the solution was a maximum or minimum of the energy landscape as a function of the stretching parameter . Here we would like to do a similar analysis to study whether the energy increases or decreases as we let where is a small translation. Again we will posit that there is a translational instability if the energy decreases as departs from zero. Let us again consider NLSE plus real potential

 iψt+ψxx+g|ψ|2κψ−V(x)ψ=0, (39)

where . We want to see how the energy of the system changes under the translation with the normalization fixed by the requirement that the mass is preserved. Clearly choosing

 ψa=ψ(x+a) (40)

preserves the mass for the wave function of the exact solution whose dependence is displayed by letting

 ψ(x)=A(b,γ)\sechγ(x)eiγ2t (41)

It is also clear that both and remain unchanged under . Only is not translationally invariant. We define

 H3(a,γ) =∫\ddxV(x)|ψa|2 (42) =λ~b2A20∫dx\sech2(x)\sech2γ(x+a) =λ~b2A20∫dy\sech2(y−a)\sech2γ(y).

We need to ensure that as a function of the energy is stationary. Since both and are independent of we only need to consider:

 \pdvH3(a,λ)a∣∣a=0 (43)

at and then calculate the second derivative at . Clearly

 \pdvHβa=\pdvH3a(a)= (44) −2λ~b2A20∫\ddy\sech2(y−a)tanh(y−a)\sech2γ(y),

which is indeed zero when evaluated at , being an odd function of . The second derivative gives

 \pdv[2]Ha(a)=\pdv[2]H3aa= (45) 2λ~b2A20∫\ddy\sech2γ(y)[3\sech4(y−a)−2\sech2(y−a)].

Evaluating this at we find:

 \pdv[2]Haa∣∣a=0=2λ~b2A20√πγΓ[γ+1]Γ[γ+5/2]. (46)

Thus we indeed find that the solitary wave has translational instability if , i.e. for the repulsive potential, while it is stable in case , i.e. attractive potential. Note that the answer does not depend on the sign of . We show in Fig. b the three-dimensional landscape for for a stretching and displacement shift, . One can see stability for an attractive potential (case I) but a saddle point for a repulsive potential (case II).

## V Linear Stability Analysis

The traditional way to study stability under small perturbations is to perform a linear stability analysis which we now present. The results obtained agree with the simpler analysis using Derrick’s theorem and looking at the effects of translation on the energy landscape. Starting with the NLSE equation in an external potential:

 {i∂t+∂2x+g\absψ(x,t)2κ−V(x)}ψ(x,t)=0, (47)

which has solitary wave solutions , with . Here satisfies

 {ω+∂2x+g\absϕω(x)2κ−V(x)}ϕω(x)=0. (48)

For , one has the explicit expression:

 ϕ−γ2(x)=A0(~b,γ)\sechγ(x). (49)

We consider perturbations in the form with and linearize the equation (47) with respect to . The linearized equation is of the form

 {ω+i∂t+∂2x+V(x)+g\absϕω(x)2κ}r(x,t) (50) +2κg\absϕω(x)2κRr(x,t)=0.

Because of the term , Eq. (50) is not a -linear operator. For computational convenience we separate the real and imaginary parts of and define

 R(x,t)=(p(x,t)q(x,t))=(Rr(x,t)Ir(x,t)). (51)

The equation (50) can then be written as

 ∂tR=JLR (52)

where

 J=(01−10)\qcL=(L+00L−) (53)

 L−(ω) =−∂2x+V(x)−ω−g\absϕω(x)2κ, (54a) L+(ω) =L−(ω)−2κg\absϕω(x)2κ. (54b)

From (48) and its derivative with respect to , we find

 L−(ω)ϕω(x)=0\qcL+(ω)∂ωϕω(x)=ϕω(x). (55)

To explore the linear stability we consider eigenvalues for the operator . Since

 (JL)2=(−L−(ω)L+(ω)00−L+(ω)L−(ω)) (56)

and the nonzero eigenvalues of and coincide, we can consider eigenvalues of the operator instead. If has eigenvalues with positive real part, so does and the solitary wave solutions are linearly unstable; otherwise, only has purely imaginary eigenvalues and the solitary wave solutions are spectrally stable.

For case (II) , we have . The amplitude is in and , so is an eigenfunction of , corresponding to zero eigenvalue. Since is positive, is nonnegative and the kernel is span by Proposition 2.8 in Weinstein (1985). Due to (55) is an eigenfunction of , corresponding to zero eigenvalue. We will show that the smallest eigenvalue of is negative. According to Ref. Sulem, Catherine and Sulem, Pierre-Louis (1999) the smallest eigenvalue of is given by,

 minσd(L−L+)=min{\expvalu,L+u\expvalu,L−1−u,u∈ker(L−)⊥}. (57)

is positive definite in , so the sign of smallest eigenvalue is decided by that of . Since is an even function, is an odd function in . Due to the identity

 ∂x(L−ϕω(x))=L+∂xϕω(x)+V′(x)ϕω(x)=0, (58)

one has that . Hence we have

 \expval∂xϕω,L+∂xϕω=\expval∂xϕω,−V′(x)ϕω (59) =12\expvalϕω,V′′(x)ϕω≈−0.455<0.

It follows that can be negative and thus has at least one positive eigenvalue. We conclude that the solitary wave solutions are linearly unstable.

For case (III) , we know that is nonnegative as in case (II). Since is positive, is positive as well. It implies that has only negative eigenvalues and the solitary wave solutions are spectrally stable.

## Vi Four parameter time-dependent trial wave function

In order to study the dynamics in a collective coordinate approximation as well as determine the frequency of small oscillations (or the intitial growth of instabilities) we will now consider a four-parameter trial wave function of the form:

 ~ψ(x,t)=A(t)\sechγ[β(t)y(x,t)]ei~ϕ(x,t), (60)

where

 ~ϕ(x,t)=−θ(t)+p(t)y(x,t)+Λ(t)y2(x,t). (61)

Here we have put . The parameter is related to the canonical conjugate variable to the average value of . It arises naturally in the Hartree-Fock approximation to the dynamics of the Schrödinger equation Cooper et al. (1986).

It will be useful to define a reciprocal width parameter and use this parameter as a generalized coordinate. Conservation of probability gives the “mass” equation,

 M=∫\ddxρ(x,t)=G(t)A2(t)c1[γ], (62)

where is given in Eq. (10). In order to maintain probability conservation, we want to keep constant. That is, we put

 A2(t)=MG(t)c1[γ] (63)

so and are not independent variables. The phase does not enter into the Hamiltonian, and we will ignore it in what follows. The trial wave function we will assume is:

 ~ψ(x,t)=√MGc1[γ]\sechγ(y/G)ei[py+Λy2]. (64)

The four variational parameters are labeled by

 Qi(t)=\qtyq(t),p(t),G(t),Λ(t). (65)

Taking the appropriate derivatives we find that

 L0=i2∫\ddx[~ψ∗~ψt−~ψ∗t~ψ]=πi(Q)˙Qi, (66)

where

 πq =pM, (67a) πp =0 (67b) πβ =0, (67c) πΛ =−MG2c2[γ]/c1[γ], (67d)

with given in Eq. (35). The only non-zero derivatives of the are

 ∂pπq=M\qc∂GπΛ=−2MGc2[γ]/c1[γ], (68)

so the symplectic tensor is

 fij(Q)=M⎛⎜ ⎜ ⎜⎝0−1001000000−C00C0⎞⎟ ⎟ ⎟⎠\qcC=2Gc2[γ]c1[γ], (69)

and the inverse is

 fij(Q)=1MC⎛⎜ ⎜ ⎜⎝0C00−C000000100−1⎞⎟ ⎟ ⎟⎠. (70)

We also find that for the four-parameter trial wave function

 H=H1+H2+H3, (71)

where

 H1 =∫\ddx|~ψx|2 (72a) =Mp2+MG2γ2c1[γ+1]c1[γ]+4MG2Λ2c2[γ]c1[γ], H2 =−gκ+1∫\ddx|~ψ|2κ+2 (72b) =−gMγγ+1\qty(MGc1[γ])1/γc1[γ+1]c1[γ] H3 =∫\ddxV(x)|~ψ|2=−λ~b2Mf1[G,q,γ]c1[γ], (72c)

with

 f1[G,q,γ]=∫\ddz\sech2γ(z)\sech2(Gz+q). (73)

The Hamiltonian then becomes

 H(Q)M=p2+γ2G2c1[γ+1]c1[γ]+4G2Λ2c2[γ]c1[γ] (74) −gγγ+1\qty(MGc1[γ])1/γc1[γ+1]c1[γ]−λ~b2f1[G,q,γ]c1[γ].

Derivatives of the Hamiltonian, , are given by

 vqM =2λ~b2f2[G,q,γ]c1[γ], (75a) vpM =2p, (75b) vGM =−γG3c1[γ+1]c1[γ]+8GΛ2c2[γ]