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

TL;DR

This paper analyzes the stability of exact solutions to the nonlinear Schrödinger equation with Pöschl-Teller potential, revealing conditions for instability and the effects of perturbations through analytic and numerical methods.

Contribution

It provides new analytic stability criteria for exact solutions of the NLSE with Pöschl-Teller potential, considering both attractive and repulsive cases, using energy landscape and variational methods.

Findings

Repulsive potential causes translational instability leading to wave splitting.

Attractive potential solutions are more stable under small perturbations.

Numerical simulations confirm analytic stability predictions.

Abstract

We discuss the stability properties of the solutions of the general nonlinear \Schrodinger\ equation (NLSE) in 1+1 dimensions in an external potential derivable from a parity-time ($\PT$) symmetric superpotential $W(x)$ that we considered earlier \cite{PhysRevE.92.042901}. In particular we consider the nonlinear partial differential equation $ \{ i \, \partial_t + \partial_x^2 - V(x) + g | \psi(x,t) |^{2\kappa} \} \, \psi(x,t) = 0 \>, $ for arbitrary nonlinearity parameter $\kappa$, where $g= \pm1$ and $V$ is the well known P{\"o}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 $\kappa$. For the repulsive potential (and $g=+1$) 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 $g=+1$), changing the initial conditions by 1 \% modifies the domain of stability only slightly. For the case of the attractive potential and negative $g$ perturbed solutions merely oscillate with the oscillation frequencies predicted by the variational approximation.

Figures

Equations

$$\{\mathrm{i}{\partial }_t+{\partial }_x^2+g\text{absolutevalue}{\psi (x,t)}^{2\kappa }-V(x)\}\psi (x,t)=0,$$

$$V(x)=-\lambda {\tilde{b}}^2{\text{sech}}^2(x),{\tilde{b}}^2=b^2-1/4,$$

$$W(x)=\text{quantity}(m-1/2)\tanh x-ib\text{sech}x,$$

$${\psi }_0(x,t)$$

$$A_0^{2/\gamma }(\tilde{b},\gamma )$$

$${\tilde{b}}_{\gamma }^2\equiv \gamma (\gamma +1)\ge {\tilde{b}}^2\ge 0.$$

$${\psi }_0(x,t)$$

$$A_0^{2/\gamma }(\tilde{b},\gamma )$$

$${\psi }_0(x,t)$$

$$A_0^{2/\gamma }(\tilde{b},\gamma )$$

$$g{A}_0^{2/\gamma }(\tilde{b},\gamma )=\gamma (\gamma +1)-\lambda {\tilde{b}}^2,$$

$$M_0(\tilde{b},\gamma )=\int \text{differential}x\text{absolutevalue}{{\psi }_0(x,t)}^2=A_0^2(\tilde{b},\gamma )c_1[\gamma ],$$

$$c_1[\gamma ]=\int \text{differential}z{\text{sech}}^{2\gamma }(z)=\frac{\sqrt{\pi }\Gamma [\gamma ]}{\Gamma [\gamma +1/2]}.$$

$$\Gamma [\psi ,{\psi }^{\ast }]=\int \text{differential}tL[\psi ,{\psi }^{\ast }]$$

$$L[\psi ,{\psi }^{\ast }]$$

$$H[\psi ,{\psi }^{\ast }]$$

$$\frac{\delta \Gamma }{\delta {\psi }^{\ast }}=\frac{\delta \Gamma }{\delta \psi }=0.$$

$$\tilde{\psi }[x,Q(t)],Q(t)=\{Q^1(t),Q^2(t),\dots ,Q^{2n}(t)\}.$$

$$L[Q,\dot{Q}]={\pi }_i(Q){\dot{Q}}^i-H[Q],$$

$${\pi }_i(Q)=\frac{\mathrm{i}}{2}\int \text{differential}x\{{\tilde{\psi }}^{\ast }[{\partial }_i\tilde{\psi }]-[{\partial }_i{\tilde{\psi }}^{\ast }]\tilde{\psi }\},$$

$$H(Q)=\!\int\!\differential{x}\bigl{[}\,|\partial_{x}\tilde{\psi}|^{2}-\frac{g\,|\tilde{\psi}|^{2\kappa+2}}{\kappa+1}+V(x)\,|\tilde{\psi}|^{2}\,\bigr{]}\>.$$

$$\derivative{t}\Bigl{(}\partialderivative{L}{\dot{Q}^{i}}\Bigr{)}-\partialderivative{L}{Q^{i}}=0\>.$$

$$f_{ij}(Q){\dot{Q}}^j=v_i(Q),$$

$$f_{ij}(Q)={\partial }_i{\pi }_j(Q)-{\partial }_j{\pi }_i(Q)$$

$$f^{ij}(Q)f_{jk}(Q)={\delta }_k^i,$$

$${\dot{Q}}^i=f^{ij}(Q)v_j(Q).$$

$$\text{derivative}H(Q)t={\dot{Q}}^i{v}_i(Q)=f^{ij}(Q)v_j(Q)v_i(Q)=0,$$

$$\{A(Q),B(Q)\}=({\partial }_{i}A(Q))f^{ij}(Q)({\partial }_{j}B(Q)).$$

$$\{Q^i,Q^j\}=f^{ij}(Q).$$

$${\psi }_{\beta }(x,t)=A(\tilde{b},\beta ,\gamma ){\text{sech}}^{\gamma }(\beta x),$$

$$M[\tilde{b},\beta ,\gamma ]$$

$$=M_0[\tilde{b},\gamma ]=A_0^2(\tilde{b},\gamma )c_1[\gamma ].$$

$$H(\beta ,\gamma )=H_1(\beta ,\gamma )+H_2(\beta ,\gamma )+H_3(\beta ,\gamma ),$$

Full text

††preprint: NLSE-Stability-v10.tex, , \now EDT

I Introduction

The study of open systems with balanced loss and gain, typically defined by Parity-Time ($\mathcal{PT}$) 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 $\mathcal{PT}$-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, $\mathcal{PT}$-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 $\mathcal{PT}$-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, $\mathcal{PT}$-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 $\mathcal{PT}$-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 $g$ 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 $g=\pm 1$ and the sign of the potential $\lambda=\pm 1$ we have found different classes of exact soltuions when the NLSE is in the presence of a Pöschl-Teller potential centered at $x=0$. Schrödinger’s equation for these cases is given by:

[TABLE]

where

[TABLE]

with $\tilde{b}^{2}>0$ and $\kappa>0$. Since $g$ can be scaled out of the equation by letting $\psi(x,t)\mapsto g^{-1/(2\kappa)}\,\psi(x,t)$, we can restrict ourselves to $g=\pm 1$ in what follows. The signs here are chosen such that the nonlinear term is attractive for $g=+1$ and repulsive for $g=-1$ and the external Pöschl-Teller potential is attractive for $\lambda=+1$ and repulsive for $\lambda=-1$. This potential is a special case of potentials obtainable from the complex $\mathcal{PT}$-symmetric SUSY superpotential

[TABLE]

with $m=1$, which gives rise to $\mathcal{PT}$-symmetric partner potentials $V_{\pm}=W^{2}\pm W^{\prime}$. Our real $V(x)$ corresponds to $V_{+}$. There are several cases of Eq. (1) which have exact solutions. These are

1. (I)

Attractive nonlinear term and attractive potential: $g=+1$, $\lambda=+1$. In this case, the exact solution is given by

[TABLE]

where $\gamma=1/\kappa$. In this case,

[TABLE]

We studied this case in a previous paper Cooper et al. (2017), where we found that all solitary waves for $\kappa<2$ and $0<\tilde{b}^{2}<\tilde{b}^{2}_{\gamma}$ are stable, as for the case of solitary waves in the NLSE ($\tilde{b}^{2}=0$). However, we also found a new region above $\kappa=2$ where these solutions are stable. 2. (II)

Attractive nonlinear term and repulsive potential: $g=+1$, $\lambda=-1$. For this case, the exact solution is given by

[TABLE]

In this case, we only require $\tilde{b}^{2}\geq 0$. This solution goes over to a particular moving solitary wave solution of the NLSE when $\tilde{b}\rightarrow 0$. Since the solutions of the NLSE are stable to deformations of the width for all $\kappa<2$ we expect (and we will find) that in that regime there will be a critical value of $\tilde{b}$ above which the solution will be unstable to width deformations. We expect and we find that again the solutions are always unstable for $\kappa>2$. What we will also find is that for all values of $\kappa$ these solutions are unstable to a slight translation, even if induced by numerical noise. 3. (III)

Repulsive nonlinear term and attractive potential: $g=-1$, $\lambda=+1$. For this case, the exact solution is given by

[TABLE]

In this case, we require $\tilde{b}^{2}\geq\tilde{b}^{2}_{\gamma}$. For this choice of $g$ 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

[TABLE]

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

[TABLE]

where

[TABLE]

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.111P. A. M. Dirac, Appendix to the Russian edition of The Principles of Wave Mechanics, as cited by Ia. I. Frenkel, Wave Mechanics, Advanced General Theory (Clarendon Press, Oxford, 1934), pp. 253, 436. Pattanayak and Schieve Pattanayak and Schieve (1994) point out that the reference often quoted, P. A. M. Dirac, Proc. Cambridge Philos. Soc. 26, 376 (1930), does not contain this equation. 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

[TABLE]

where

[TABLE]

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

[TABLE]

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:

[TABLE]

Introducing the notation $\partial_{i}\equiv\partial/\partial Q^{i}$, the Lagrangian (12a) is given by

[TABLE]

where

[TABLE]

and $H(Q)$ is given by

[TABLE]

The Euler-Lagrange equations now become

[TABLE]

From (14) this gives

[TABLE]

where we have set $v_{i}(Q)\equiv\partial_{i}H(Q)$, and where

[TABLE]

is an antisymmetric $2n\times 2n$ symplectic matrix. If $\det{f(Q)}\neq 0$, we can define an inverse as the contra-variant matrix with upper indices,

[TABLE]

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

[TABLE]

Conservation of energy is expressed as

[TABLE]

since $f^{ij}(Q)$ is an antisymmetric tensor. Poisson brackets are defined using this antisymmetric tensor. If $A(Q)$ and $B(Q)$ are functions of $Q$, Poisson brackets are defined by

[TABLE]

In particular,

[TABLE]

This definition satisfies Jacobi’s identity. That is, what we have shown here is that the $2n$ quantities $Q^{i}$ 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, $x\mapsto\beta x$ with fixed normalization, the energy of the system is lowered. The stretched wave function for Derrick’s theorem is given by

[TABLE]

with the normalization fixed by the requirement,

[TABLE]

So $A^{2}(\tilde{b},\beta,\gamma)=\beta A_{0}^{2}(\tilde{b},\gamma)$. Evaluation of the Hamiltonian (12b) with Derrick’s wave function gives:

[TABLE]

where

[TABLE]

where

[TABLE]

So then $H(\beta,\gamma)/M_{0}[\tilde{b},\gamma]$ is given by

[TABLE]

and is independent of $g$. $h(\tilde{b},\beta,\gamma)$ is stationary when $\beta=1$. We have

[TABLE]

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

[TABLE]

so that at $\beta=1$,

[TABLE]

for all values of $\gamma$ and $\tilde{b}^{2}$. The sign of the second derivative of $h(\tilde{b},\beta,\gamma)$ with respect to $\beta$ at $\beta=1$ determines whether the solution is unstable to small changes in the width. If $\partial^{2}h(\tilde{b},\beta,\gamma)/\partial\beta^{2}$ evaluated at $\beta=1$ is negative the solution is unstable. We find

[TABLE]

From Appendix A, we have

[TABLE]

where

[TABLE]

Inserting this into (34) and evaluating it at $\beta=1$ gives:

[TABLE]

which is independent of $g$.

The critical value of $\tilde{b}^{2}$, when (36) vanishes, is

[TABLE]

which is independent of $g$. One can easily check that

[TABLE]

for all $\gamma$. In Appendix A, we give an alternative form for Eq. (37), which is in agreement with Ref. [Cooper et al., 2017]. In Fig. 1, we plot $\tilde{b}^{2}_{\text{crit}}$ for the three cases. Referring to the figure, according to Derrick’s theorem,

1. (I)

For case I with $g=+1$, $\lambda=+1$, and $0\leq\tilde{b}^{2}\leq\tilde{b}^{2}_{\gamma}$, we see from Eq. (36) that $\partial^{2}h(\tilde{b},\beta,\gamma)/\partial\beta^{2}\geq 0$ for all $\kappa<2$, so Derrick’s theorem predicts that solutions are width stable for $\kappa<2$. For $\kappa>2$, there is another region for $\tilde{b}^{2}_{\text{crit}}<\tilde{b}^{2}<\tilde{b}^{2}_{\gamma}$ where width stable solutions are also possible. 2. (II)

For case II with $g=+1$, $\lambda=-1$, and $0\leq\tilde{b}^{2}\leq\tilde{b}^{2}_{\text{crit}}$, width stable solutions are possible for $\kappa<2$, but we will find that this region is unstable to translation perturbations. 3. (III)

For case III with $g=-1$, $\lambda=+1$, and $\tilde{b}^{2}_{\gamma}<\tilde{b}^{2}$ width stable solutions are possible for all $\kappa$.

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 $\beta$, 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 $\beta$. Here we would like to do a similar analysis to study whether the energy increases or decreases as we let $x\rightarrow x+a$ where $a$ is a small translation. Again we will posit that there is a translational instability if the energy $H[a]$ decreases as $a$ departs from zero. Let us again consider NLSE plus real potential

[TABLE]

where $V(x)=-\lambda\tilde{b}^{2}\sech^{2}(x)$. We want to see how the energy of the system changes under the translation $x\rightarrow x+a$ with the normalization fixed by the requirement that the mass $M$ is preserved. Clearly choosing

[TABLE]

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

[TABLE]

It is also clear that both $H_{1}$ and $H_{2}$ remain unchanged under $x\rightarrow x+a$. Only $H_{3}$ is not translationally invariant. We define

[TABLE]

We need to ensure that as a function of $a$ the energy is stationary. Since both $H_{1}$ and $H_{2}$ are independent of $a$ we only need to consider:

[TABLE]

at $\beta=1$ and then calculate the second derivative at $\beta=1$. Clearly

[TABLE]

which is indeed zero when evaluated at $a=0$, being an odd function of $y$. The second derivative gives

[TABLE]

Evaluating this at $a=0$ we find:

[TABLE]

Thus we indeed find that the solitary wave has translational instability if $\lambda<0$, i.e. for the repulsive potential, while it is stable in case $\lambda>0$, i.e. attractive potential. Note that the answer does not depend on the sign of $g$. We show in Fig. 2 the three-dimensional landscape for $H(\beta,a)$ for a stretching and displacement shift, $x\rightarrow\beta x+a$. 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:

[TABLE]

which has solitary wave solutions $\psi(x,t)=\phi_{\omega}(x)\,{\rm e}^{-{\rm i}\omega t}$, with $\phi_{\omega}(x)\in\mathbb{R}$. Here $\phi_{\omega}(x)$ satisfies

[TABLE]

For $\omega=-\gamma^{2}$, one has the explicit expression:

[TABLE]

We consider perturbations in the form $\psi(x,t)=(\phi_{\omega}(x)+r(x,t)){\rm e}^{-{\rm i}\omega t}$ with $r(x,t)\in\mathbb{C}$ and linearize the equation (47) with respect to $r(x,t)$. The linearized equation is of the form

[TABLE]

Because of the term $\Re{r(x,t)}$, Eq. (50) is not a $\mathbb{C}$-linear operator. For computational convenience we separate the real and imaginary parts of $r(x,t)$ and define

[TABLE]

The equation (50) can then be written as

[TABLE]

where

[TABLE]

with self-adjoint operators

[TABLE]

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

[TABLE]

To explore the linear stability we consider eigenvalues for the operator $\mathbf{JL}$. Since

[TABLE]

and the nonzero eigenvalues of $(-L_{-}L_{+})$ and $(-L_{+}L_{-})$ coincide, we can consider eigenvalues of the operator $(-L_{-}L_{+})$ instead. If $(-L_{-}L_{+})$ has eigenvalues with positive real part, so does $\mathbf{JL}$ and the solitary wave solutions are linearly unstable; otherwise, $\mathbf{JL}$ only has purely imaginary eigenvalues and the solitary wave solutions are spectrally stable.

For case (II) $g=+1$, $\lambda=-1$ we have $\omega=-\gamma^{2}$. The amplitude $\phi_{\omega}(x)$ is in $L^{2}$ and $L_{-}(\omega)\phi_{\omega}=0$, so $\phi_{\omega}$ is an eigenfunction of $L_{-}(\omega)$, corresponding to zero eigenvalue. Since $\phi_{\omega}$ is positive, $L_{-}(\omega)$ is nonnegative and the kernel $\ker(L_{-}({\omega}))$ is span$\{\phi_{\omega}\}$ by Proposition 2.8 in Weinstein (1985). Due to (55) $\partial_{\omega}\phi_{\omega}$ is an eigenfunction of $L_{-}L_{+}$, corresponding to zero eigenvalue. We will show that the smallest eigenvalue of $L_{-}L_{+}$ is negative. According to Ref. Sulem, Catherine and Sulem, Pierre-Louis (1999) the smallest eigenvalue of $L_{-}L_{+}$ is given by,

[TABLE]

$L_{-}$ is positive definite in $\ker(L_{-})^{\perp}$, so the sign of smallest eigenvalue is decided by that of $\langle u,L_{+}u\rangle$. Since $\phi_{\omega}(x)$ is an even function, $\partial_{x}\phi_{\omega}(x)$ is an odd function in $\ker(L_{-})^{\perp}$. Due to the identity

[TABLE]

one has that $L_{+}\partial_{x}\phi_{\omega}(x)=-V^{\prime}(x)\phi_{\omega}(x)$. Hence we have

[TABLE]

It follows that $\expectationvalue{u,L_{+}u}$ can be negative and thus $(-L_{-}L_{+})$ has at least one positive eigenvalue. We conclude that the solitary wave solutions are linearly unstable.

For case (III) $g=-1$, $\lambda=+1$ we know that $L_{-}$ is nonnegative as in case (II). Since $-2\kappa g\absolutevalue{\phi_{\omega}(x)}^{2\kappa}$ is positive, $L_{+}$ is positive as well. It implies that $-L_{-}L_{+}$ 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:

[TABLE]

where

[TABLE]

Here we have put $y(x,t)=x-q(t)$. The parameter $\Lambda$ is related to the canonical conjugate variable to the average value of $y^{2}$. 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 $G(t)=1/\beta(t)$ and use this parameter as a generalized coordinate. Conservation of probability gives the “mass” equation,

[TABLE]

where $c_{1}[\gamma]$ is given in Eq. (10). In order to maintain probability conservation, we want to keep $M$ constant. That is, we put

[TABLE]

so $A(t)$ and $G(t)$ are not independent variables. The phase $\theta(t)$ does not enter into the Hamiltonian, and we will ignore it in what follows. The trial wave function we will assume is:

[TABLE]

The four variational parameters are labeled by

[TABLE]

Taking the appropriate derivatives we find that

[TABLE]

where

[TABLE]

with $c_{2}[\gamma]$ given in Eq. (35). The only non-zero derivatives of the $\pi_{i}$ are

[TABLE]

so the symplectic tensor is

[TABLE]

and the inverse is

[TABLE]

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

[TABLE]

where

[TABLE]

with

[TABLE]

The Hamiltonian then becomes

[TABLE]

Derivatives of the Hamiltonian, $v_{i}=\partial_{i}H(Q)$, are given by

[TABLE]

where

[TABLE]

From Eq. (21) and using (70), we find

[TABLE]

Recall that from (63),

[TABLE]

and since $G(0)=1$ at $t=0$, we can use Eq. (7) to find:

[TABLE]

So the dynamic equations (77) reduce to:

[TABLE]

and are independent of $g$. We need to remember that for case III when $g=-1$ and $\lambda=+1$ there is no solution unless $\tilde{b}^{2}>\gamma(\gamma+1)$.

VI.1 Small oscillations

We require that the variational wave function starts out so that it agrees with the exact solutions at $t=0$,

[TABLE]

with $A_{0}$ fixed by (79). This means that we want to choose

[TABLE]

Setting

[TABLE]

to first order, we have

[TABLE]

where

[TABLE]

where in the last term we have integrated by parts, and

[TABLE]

As a consistency check, at $t=0$ Eqs. (80) must all give values consistent with the initial conditions. It is obvious that the first three equations are consistent, namely they all evaluate to zero. At $t=0$, Eq. (80d) becomes

[TABLE]

as required. Substituting Eqs. (VI.1) into (77) gives

[TABLE]

So to first order, the $(q,p)$ and $(G,\Lambda)$ modes uncouple and reduce to equations of the form:

[TABLE]

where the $(q,p)$ mode frequency is given by

[TABLE]

independent of $g$. For $\lambda=+1$, translational motion is stable for all $\gamma$, and unstable for $\lambda=-1$. This is easily explained by the fact that as long as the solitary wave is near $q=0$, for $\lambda=+1$ it sees an attractive force that brings it back to the origin. For the opposite sign, $\lambda=-1$, it sees a repulsive force that moves it from the origin, assuming it maintains its shape. We illustrate this behavior in Sec. VII.2 by a numerical solution of Schrödinger’s equation using a split-operator method.

The $(G,\Lambda)$ mode frequency is given by

[TABLE]

where

[TABLE]

and the critical value of $\tilde{b}^{2}$ is then given by

[TABLE]

which is the same result for $\tilde{b}_{\text{crit}}^{2}$ that we found in Eq. (37) using Derrick’s theorem. Here in addition, we find a value for the $(G,\Lambda)$ mode frequency $\omega_{G}^{2}$ in Eq. (91). Plots of $\tilde{b}_{\text{crit}}^{2}$, which are the same as Derrick’s theorem, are shown in Figs. 1.

VI.2 Dynamics of the Collective coordinates

In this section we plot representative time evolutions of the collective coordinates of the four-parameter variational calculation given in Eqs. (80) for the three cases. The parameters $(\,q,p,G,\Lambda\,)$ can be related to the expectation vales of $\langle\,x,p_{\text{op}},x^{2},xp_{\text{op}}\,\rangle$ where $p_{\text{op}}=-i\partial/\partial x$. Here $\langle x\rangle/M=q(t)$, etc. One can show that using the equations obeyed by these four collective coordinates, that the evolution equations for these expectation values are exactly satisfied. We expect, and find that these time evolutions qualitatively agree with the numerically calculated values of these collective coordinates especially when the actual form of the numerically determined wave function is preserved. In case III, although we predict the translational instability quite well, we did not anticipate that the form of the wave function would bifurcate. Of course a posteriori one could assume a two humped variational wave function with more parameters to actually capture better the time evolution of the initial exact solution. What we find, is that when we are in the oscillatory regime of either $q(t)$ or $G(t)$, the oscillation period determined from the small oscillation equations is in good agreement with what is found from the dynamical evolution of $q(t)$ and $G(t)$ from their evolution equations.

VI.2.1 $g=+1$ and $\lambda=+1$ (Case I)

This is the case that we studied in our previous paper Cooper et al. (2017). However, in that paper we did not consider oscillations in the spatial direction, nor did we numerically solve the NLSE to compare with our analytic results. For this case in the $G$ unstable region the behavior in $q$ is oscillatory in our approximation. We will choose $\kappa=5/2$ and two values of $\tilde{b}^{2}=1/10$ which is in the unstable regime and $\tilde{b}^{2}=1/5$ which is in the stable regime to display the two types of behavior for the parameters $q(t)$ and $G(t)$ as a function of time. For the blowup case, the period for $q(t)$ is oscillatory with a period of $T=32.5$, which agrees with the numerical results in Fig. 3(a). However $G(t)$ blows up, $G(t)\rightarrow 0$, as shown in Fig. 3(b).

For $\tilde{b}^{2}=1/5$ one is in the oscillatory regime for $G(t)$ and $q(t)$ and we get the results shown in Figs. 4. For this case the small oscillation equation predicts that the period of $q(t)$ is $T=23$, and the period of $G(t)$ is $T=37$.

VI.2.2 $g=+1$ and $\lambda=-1$ (Case II)

The case where we have an exact solution for a repulsive potential leads to the most unexpected behavior, as we will show later in our numerical simulations. The first interesting thing is that although Derrick’s theorem shows that the answer is stable to changing the width when $\kappa<2$ and $\tilde{b}^{2}$ is below $\tilde{b}^{2}_{\text{c}}$ in Fig. 1(b), we find that if we shift the position by a small amount, because of the repulsive potential, the solitary wave is pushed out of the region of the potential and then oscillates about a potential free solution in this approximation. This result is suggested by the four-component variational calculation, and confirmed by a numerical calculation shown in Fig. 9(b). Choosing $\kappa=1$ and $\tilde{b}^{2}=1$, which is in the regime which is stable to width changes, the solutions of the dynamic equations (78) give the results shown in Fig. 5.

VI.2.3 $g=-1$ and $\lambda=+1$ (Case III)

For this case, all allowed solutions $\tilde{b}^{2}>(\kappa+1)/\kappa^{2}$ should be stable to small changes in both $G$ and $q(t)$. This oscillatory behavior for the case $\kappa=3$, $\tilde{b}^{2}=1/2$, $q[0]=0.01$ is shown in Figs. 7.

VII Numerical Study of Stability

VII.1 Domains of stability

In order to study the stability of the soliton solutions, the actual numerical simulation of the soliton evolution has been performed. For that purpose, we have numerically solved Eq. (47) with the initial conditions described in Section I.1 using a Crank-Nicolson scheme Press et al. (1992).

The complex soliton in the spatial domain was represented on a regular grid with mesh size $\Delta x=5\times 10^{-3}$, and free boundary conditions were imposed.

In order to study the stability regimes, we calculate the normalized correlation between the initial intensity profile (at $t=0$), and the intensity profile for $t>0$, i.e.

[TABLE]

Notice that $0\leq\text{corr}(0,t)=C_{t}<\infty$ and its value $C_{t}$ can be interpreted as follows

[TABLE]

for any $t>0$. Theoretically, it means that the evolution of the soliton solutions should be checked up to $t\to\infty$. However, in practice, the inherent numerical noise of the simulations randomly perturbs the soliton shape during evolution, so finite evolution times are enough to determine the stability of soliton solutions. In this regard, we have found that the evolution of solitons up to $5\times 10^{2}$ time units with step size $\Delta t=3\times 10^{-2}$ is enough to study their stability. First let us return to the problem we studied earlier (Case I) where $\lambda=g=1$. For that case three methods predicted the same width stability region, namely Derrick’s theorem, setting $\omega_{G}^{2}=0$ and the V-K criterion. In this case, if we consider the domain of numerical stability using the exact initial conditions, we agree with the results of Derrick’s theorem.

Next look at the case when our variational method predicts stability, namely for the attractive potential with $g=-1$. The stability should occur as long as there is a solution namely $\tilde{b}^{2}>(\kappa+1)/\kappa^{2}$ as shown in Fig. (1a). The numerical solutions lead to the same conclusions as shown in Fig. 8(a). Figures. 8(a) and 8(b) show the distribution of the stability region, blow-up regime, and translational instability for attractive and repulsive potentials, respectively.

In Fig. 8(b) we show the regions of translational and width instabilities for Case II.

VII.2 Effect of Translational Instability

This is illustrated for parameters in region II in Fig. 9 by numerical solutions of Schrödinger’s equation using a split-operator method. In Fig. 9(a), the initial conditions are such that the solution “rolls” off the top of the potential in both directions and bifurcates, reminiscent of the quantum roll problem discussed in connection with the inflationary universe Cooper et al. (1986). Although classically if we are at a maximum or saddle point we expect a particle to go to the right or to the left, since the underlying theory is essentially quantum mechanical, some of the mass density can go to the right or to the left. If we start at the top of the saddle, and the repulsive potential is strong ($\tilde{b}^{2}\sim 1$) then the solution bifurcates. In Fig. 9(b), we give a slight kick away from the origin by setting $q(0)=0.01$, which leads to translational motion away from the origin. In this region of parameter space, the results are strongly dependent on the initial conditions.

VII.3 Effects of initial conditions

A small perturbation of the initial conditions can also lead to instabilities. These effects smear out the instability regions predicted by Derrick’s theorem as shown in Fig. 1. As an illustration of such effects, we show in Fig. 10 the modification of region I for $\kappa>2$ as a result of a 1% change in $G(0)$ (using $G(0)=1.01$) instead of the exact solution value of $G(0)=1.00$. The unstable regime found by a numerical solution of Schrödinger’s equation with this initial condition is somewhat bigger than that found from Derrick’s theorem (see Fig. 10). Similar effects are found for other initial conditions.

VIII Conclusions

In this paper we studied the stability of two new exact classes of solutions of the NLSE in a real Pöschl-Teller potential which is the SUSY partner of a complex $\mathcal{PT}$-symmetric potential studied previously Kevrekidis et al. (2015). Since the exact equations are derivable from an action principle, if we approximate the wave function by a set of collective coordinates we obtain a symplectic formulation of the dynamics of the collective coordinates with a conserved Hamiltonian. This reduced system is amenable to several approaches to studying the stability problem such as Derrick’s theorem, looking at the energy landscape as a function of translations, as well as a time dependent variational approach. Using these methods we have mapped the domain of stability of these exact solutions and found good agreement with what we obtain from numerical simulations. For the case of the solution with an attractive potential but $g=-1$, the solutions are stable to both small width and small position deformations for all $\kappa$ as long as $\tilde{b}^{2}>\gamma(\gamma+1)$. In that situation, the small oscillation equations of the variational approach give good agreement with numerical simulations. For the more interesting case of $g=1$ and a repulsive potential, there is a translational instability. This can be seen by looking at the energy landscape or by looking at the small oscillation equations of our variational approximation. The stability results obtained from looking at the energy landscape of the solution as a function of both width stretching and translations are in agreement with the results of a more rigorous linear stability analysis. For the case of the repulsive potential which has a translational instability we find by performing numerical simulations quite interesting results. If one starts with the exact initial solution, the solution breaks into two equal amplitude pulses going in opposite directions. By perturbing the solution in one direction, the majority of the wave goes in that direction, but still some of the wave goes in the opposite direction.

Acknowledgements.

F.C. would like to thank the Santa Fe Institute and the Center for Nonlinear Studies at Los Alamos National Laboratory for their hospitality. A.K. is grateful to Indian National Science Academy (INSA) for awarding him INSA Senior Scientist position at Savitribai Phule Pune University, Pune, India. B.M. and J.F.D. would like to thank the Santa Fe Institute for their hospitality. B.M. acknowledges support from the National Science Foundation through its employee IR/D program. The work of A.S. was supported by the U.S. Department of Energy. E.A. gratefully acknowledges support from the Fondo Nacional de Desarrollo Científico y tecnológico (FONDECYT) project No. 1141223 and from the Programa Iniciativa Científica Milenio (ICM) Grant No. 130001. The research of A.C. was carried out at the Institute for Information Transmission Problems, Russian Academy of Sciences at the expense of the Russian Foundation for Sciences (Project 14-50-00150). The work of R.L. is supoported by the Department of Mathematics at Texas A&M University.

Appendix A Useful integrals and definitions

[TABLE]

A useful result is

[TABLE]

In Section III, we defined an integral,

[TABLE]

The first derivative of $g_{1}[\beta,\gamma]$ with respect to $\beta$ evaluated at $\beta=1$ is given by

[TABLE]

where we have integrated by parts. The second derivative, evaluated at $\beta=1$ is

[TABLE]

The last term can be written as

[TABLE]

So (99) becomes

[TABLE]

Using (35), Eq. (100) becomes

[TABLE]

We use this result in Section III. A second form of the critical curve can be found using an identity for $c_{2}[\gamma+1]$, which we derive here. First note that

[TABLE]

where

[TABLE]

Using the identity,

[TABLE]

integrating it over $z$, and evaluating at $\lambda\rightarrow 1$ gives:

[TABLE]

Substituting in $I[\gamma]$ from Eq. (102) gives:

[TABLE]

Using this result, and after a bit of algebra, it is easy to show that the critical curve Eq. (37) can be written as

[TABLE]

This form is identical to Eq. (4.28) in Ref. [Cooper et al., 2017].