Response of exact solutions of the nonlinear Schrödinger equation to small perturbations in a class of complex external potentials having supersymmetry and parity-time symmetry

LA-UR-17-25332 Response of exact solutions of the nonlinear Schr¨odinger equation to small perturbations in a class of complex external potentials having supersymmetry and parity-time symmetry Fred Cooper,1, 2, ∗John F. Dawson,3, †Franz G. Mertens,4, ‡Edward Ar´evalo,5, §Niurka R. Quintero,6, 7, ¶Bogdan Mihaila,8, ∗∗ Avinash Khare,9, †† and Avadh Saxena2, ‡‡ 1The Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA 2Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA 3Department of Physics, University of New Hampshire, Durham, NH 03824, USA 4Physikalisches Institut, Universit¨at Bayreuth, D-95440 Bayreuth, Germany 5Pontifical Catholic University of Chile, Instituto de F´ısica, Santiago, Regi´on Metropolitana, Chile 6Instituto Carlos I de Fisica Te´orica y Computacional, Universidad de Granada, E-18015, Granada, Spain 7Departamento de Fisica Aplicada I, E.P.S. Universidad de Sevilla, 41011 Sevilla, Spain 8Physics Division, National Science Foundation, Arlington, VA 22230, USA 9Physics Department, Savitribai Phule Pune University, Pune 411007, India (Dated: July 7, 2017, 12:46am EDT) We discuss the effect of small perturbation on nodeless solutions of the nonlinear Schr¨odinger equation in 1+1 dimensions in an external complex potential derivable from a parity-time symmetric superpotential that was considered earlier [Phys. Rev. E 92, 042901 (2015)]. In particular we consider the nonlinear partial differential equation {i∂t+∂2 x+g|ψ(x, t)|2−V+(x)}ψ(x, t) = 0, where V+(x) = −b2−m2+ 1/4sech2(x)−2i m b sech(x) tanh(x) represents the complex potential. Here we study the perturbations as a function of band musing a variational approximation based on a dissipation functional formalism. We compare the result of this variational approach with direct numerical simulation of the equations. We find that the variational approximation works quite well at small and moderate values of the parameter bm which controls the strength of the imaginary part of the potential. We also show that the dissipation functional formalism is equivalent to the generalized traveling wave method for this type of dissipation. I. INTRODUCTION The topic of balanced loss and gain or parity-time (PT ) symmetry and its relevance for physical applications on the one hand, as well as its mathematical structure on the other, have drawn considerable attention from both the physics and the mathematics community. The original proposal of Bender and his collaborators [1–4] towards the study of such systems was made as an alternative to the postulate of Hermiticity in quantum mechanics. Keeping in perspective the formal similarity of the Schr¨odinger equation with Maxwell’s equations in the paraxial approximation, it was realized that such PT invariant systems can in fact be experimentally realized in optics [5–14]. Subsequently, these efforts motivated experiments in several other areas including PT invariant electronic circuits [15,16], mechanical circuits [17], and whispering-gallery microcavities [18]. Concurrently, the notion of supersymmetry (SUSY) originally espoused in high-energy physics has also been ∗coop[email protected] †john.da[email protected] ‡F[email protected] §earev[email protected] ¶niurk[email protected] ∗∗ [email protected] †† [email protected] ‡‡ av[email protected]v realized in optics [19,20]. The key idea is that from a given potential one can obtain a SUSY partner potential with both potentials possessing the same spectrum, except possibly for one eigenvalue [21,22]. Therefore, an interplay of SUSY with PT symmetry is expected to be quite rich and is indeed very useful in achieving transparent as well as one-way reflectionless complex optical potentials [23–27]. A previous paper [28] explored the interplay between PT symmetry, SUSY and nonlinearity. That paper derived exact solutions of the general nonlinear Schr¨odinger (NLS) equation in 1+1 dimensions when in a PT - symmetric complex potential [21,29]. In particular, they considered the nonlinear partial differential equation i ∂t+∂2 x−V±(x) + g|ψ(x, t)|2κψ(x, t)=0,(1.1) for arbitrary nonlinearity parameter κ, with V±(x) = W2 1(x)∓W0 1(x)−(m−1/2)2,(1.2) and the partner potentials arise from the superpotential W1(x) = (m−1/2) tanh x−ib sech x , (1.3) giving rise to V+(x) = −b2−m2+ 1/4sech2(x) (1.4a) −2i m b sech(x) tanh(x), V−(x) = −b2−(m−1)2+ 1/4sech2(x) (1.4b) −2i(m−1) bsech(x) tanh(x). arXiv:1707.01574v1 [nlin.PS] 5 Jul 2017 2 For m= 1, the complex potential V+(x) has the same spectrum, apart from the ground state, as the real potential V−(x) and this fact was used in the numerical study of the stability of the bound state solutions of the NLS equation in the presence of V+(x) (see Ref. [28]). In a recent complementary study [30] of this system of nonlinear Schr¨odinger equations in PT symmetric SUSY external potentials, the stability properties of the bound state solutions of NLS equation in the presence of the external real SUSY partner potential V−(x) were investigated. The stability regime of these solutions, which depended on the parameters (b, κ), was compared to the stability regime of the related solitary wave solutions to the NLS equation in the absence of the external potential. Because the NLS equation in the presence of V−(x) is a Hamiltonian dynamical system, in Ref. [30] they were able to use several variational methods to study the stability of the solutions when they undergo certain small deformations, and showed that these variational methods agreed with a linear stability analysis based on the Vakhitov-Kolokolov (V-K) stability criterion [31,32] as well as numerical simulations that have recently been performed. In Ref. [28] we determined the exact solutions of the equation for m= 1 for V+(x), which was complex. We studied numerically the stability properties of these solutions using linear stability analysis. We found some unusual results for the stability which depended on the value of b. What was found for m= 1 (and κ= 1) was that the eigenvalues of the linear stability matrix became complex for 0.56 < b < 1.37. At that time we had not yet formulated a variational approach for deriving the NLS equation in the presence of complex potentials. Recently we have developed such an approach and have applied it to the response of the solutions of the NLS equation to weak external complex periodic potentials. Using four variational parameters we were able to successfully predict the time evolution of these solitary waves when compared to direct numerical simulation of the NLS equation in the presence of these complex potentials [33]. Given this new tool we would like to return to the original problem of the stability of the exact solutions found in Ref. [28] and see how well this variational approach agrees with numerical simulations as a function of the strength of the dissipative part of the potential which is proportional to bm. In this paper we focus on the external potential V+(x) which is symmetric in b↔m. Here we will compare the numerical simulations with the results of our collective coordinate (CC) approximation. We will also look at the linear stability analysis that arises from studying the linearization of the CC ordinary differential equations (ODEs). For the case of a real external potential, studying the eigenvalues of this reduced stability analysis predicted the correct stability regime [34]. This paper is structured as follows. In Sec. II we review the non-hermitian SUSY model that we studied in Ref. [28] and add the self-interactions of the NLS equation to the linear model. In Sec. III we give some of the exact low order moment equations for this problem. In Sec. IV we introduce our collective coordinate approach, whereas in Sec. Vwe use a four parameter trial wave function that we considered in an earlier study of soliton behavior in complex periodic external potentials, and derive equations for the four CC’s. In Sec. VI we expand the number of CC’s to six and derive equations for the six CC’s. In Sec. VII we study the linear response theory of the six CC approximation. In Sec. VIII we present our numerical strategy for solving the NLS equation starting from a perturbed exact solution. In Sec. IX we compare the four and six CC approximations with direct numerical simulations. In Sec. Xwe present our main conclusions. Finally in Appendix A we provide the definitions of various integrals and in Appendix B we show that for this problem our variational approach is equivalent to the generalized traveling wave method [35]. II. NLS EQUATION IN THE PRESENCE OF A NON-HERMITIAN SUPERSYMMETRIC EXTERNAL POTENTIAL We were interested in studying the NLS equation in the presence of a complex external potential and were intrigued by the fact that as a result of PT symmetry, there existed complex potentials whose SUSY partners were real and had explicitly known spectra of bound states. This led us to study the external potential defined by the PT symmetric SUSY superpotential W1(x) given by Eq. (1.3). This superpotential gives rise to supersymmetric partner potentials given by Eqs. (1.4). For the case m= 1, V−(x) is the well known P¨oschl-Teller potential [36,37]. The relevant bound state eigenvalues assume an extremely simple form as E(−) n=−1 4[ 2b−2n−1 ]2.(2.1) Such bound state eigenvalues only exist when n < b − 1/2. We notice that for the ground state (n=0) to exist requires b > 1/2. The existence of a first excited state (n=1) requires b > 3/2. Here we consider the general V+(x) arising from the superpotential W1(x) depending on m, b as an external potential modifying the nonlinear Schr¨odinger equation. Rewriting the external potential given in Eq. (1.4a) as, V+(x) = V1(x)+iV2(x),(2.2) we have V1(x) = −(b2+m2−1/4) sech2(x),(2.3a) V2(x) = −2mb tanh(x) sech(x).(2.3b) Note this potential is invariant under the exchange of b and m. We are interested in the stability properties of 3 the exact solutions of the NLS equation in this external potential: {i ∂t+∂2 x+g|ψ(x, t)|2κ−[V1(x) + iV2(x) ] }ψ(x, t)=0. (2.4) This equation can be obtained from a generalized Euler- Lagrange equation using a dissipation functional [33], δΓ δψ∗=−δF δψ∗ t ,(2.5) where Γ =Zdtni 2Zdx[ψ∗ψt−ψψ∗ t]−Ho,(2.6a) H=Zdxn|ψx|2−g|ψ|2κ+2 κ+ 1 +V1(x)|ψ|2o,(2.6b) F=Zdt F = iZdxdt V2(x) [ ψtψ∗−ψ∗ tψ].(2.6c) Localized solutions to Eq. (2.4) exist for arbitrary values of κ, m, b. Here we use ψ0(x, t) to denote the exact solution to the NLS equation in the external potential, ψ0(x, t) = A0sech1/κ(x) ei[ Et+φ(x) ] (2.7) where φ(x) = 4bmκ κ+ 2 tan−1[tanh(x/2) ] ,(2.8) with E= 1/κ2, and gA2κ 0=[4b2κ2−(κ+ 2)2] [4m2κ2−(κ+ 2)2] 4κ2(κ+ 2)2.(2.9) We notice when mb = 0 the potential is real and that solutions exist for m2+b2−1/4<(κ+ 1)/κ2. There are two regimes where A2 0is positive and so a solution exists when m6= 0. This form of the solution reflects the fact that the potential V+is invariant under the interchange m↔b. In a previous paper [38] we studied the stability of these solutions for m= 0 (real external potential) and for arbitrary κ. In that paper, we also considered two other cases where exact solutions exist. For the case of g=−1 and attractive potential, for V2= 0, all the solutions that were allowed were stable. Solutions also exist for V26= 0 and are given by Eqs. (2.7) and (2.9) with g=−1. For g= 1 and a repulsive real potential we found the solutions for V2= 0 were translationally unstable. Solutions again exist when V26= 0 for this case. We will not discuss these solutions further here. Here we will confine ourselves to κ= 1, g = 1 and an attractive external potential V1and study the domain of applicability of the variational methods we have developed previously to the case of increasing the dissipation by allowing mto vary. In particular for the case we will concentrate on here (κ= 1) we have that gA2 0= (4b2−9)(4m2−9)/36 ,(2.10) so that when m2<9/4 we need that b2<9/4 for there to be a solution. Also if we confine ourselves to an attractive potential so that we avoid the known translational instability associated with repulsive potentials [38], then we also require b2+m2>1/4. Note that gA2 0is independent of g. For κ= 1 we have φ(x) = (4 m b/3) tan−1[ tanh(x/2) ] ,(2.11a) ∂xφ(x) = (2/3) m b sech(x).(2.11b) III. SOME GENERAL PROPERTIES OF THE NLS EQUATION IN COMPLEX POTENTIALS We are interested in solitary wave solutions that approach zero exponentially at ±∞. For these solutions we define the mass density ρ(x, t) = |ψ(x, t)|2, and the mass or norm M(t) as M(t) =Zdx ρ(x, t) =Zdx|ψ(x, t)|2.(3.1) In addition, we define the current as: j(x, t) = i[ψ(x, t)ψ∗ x(x, t)−ψ∗(x, t)ψx(x, t) ] .(3.2) Multiplying the NLS equation (2.4) by ψ∗(x, t) and subtracting the complex conjugate of the resulting equation, we obtain ∂ρ(x, t) ∂t +∂j(x, t) ∂x = 2V2(x)ρ(x, t).(3.3) Integrating over space, and assuming that j(+∞, t)− j(−∞, t) = 0, we find dM(t) dt= 2Zdx V2(x)ρ(x, t).(3.4) Note that Mis conserved when V2(x) = 0. If we instead multiply the NLS equation by ψ∗and add the complex conjugate of the resulting equation, we get i ( ψ∗ψt−ψψ∗ t) (3.5) =−2gρ2−ψ∗ψxx −ψψ∗ xx + 2V1(x)ρ , which when we integrate over space, leads to the virial theorem: i 2Zdx(ψ∗ψt−ψ∗ tψ)−Zdxh|ψx|2−g|ψ|4i(3.6) =Zdx V1(x)|ψ|2. The average position q(t) can be defined through the first moment of xas follows: M1(t) =Zdx x ρ(x, t) = q(t)M(t).(3.7) 4 Multiplying the continuity equation (3.3) by xand integrating over all space we find: dM1 dt= 2 P(t)+2Zdx x V2(x)ρ(x, t),(3.8) where the momentum P(t) = 1 2Zdx j(x, t) (3.9) =i 2Zdx[ψ∗(x, t)ψx(x, t)−ψ∗ x(x, t)ψ(x, t) ] . Here, we assumed that lim x→∞ xj(x, t)−lim x→−∞ xj(x, t)=0.(3.10) Assuming that the density is a function of y=x−q(t) and t, we find d dt[M(t)q(t) ] = 2 P(t)+2Zdy y V2(y+q(t)) ρ(y, t) + 2 q(t)Zdx V2(x)ρ(x−q(t), t). We recognize the last term as q(t) dM(t)/dt, so that we finally have: M(t)dq(t) dt= 2 P(t)+2Zdy y V2(y+q(t)) ρ(y, t).(3.11) Taking the time derivative of the momentum P(t), using the equations of motion for ψand ψ∗, and integrating by parts, we find dP(t) dt=−Zdx ρ(x, t)∂V1(x) ∂x +Zdx j(x, t)V2(x). (3.12) Here ∂V1(x) ∂x = 2 b2+m2−1/4tanh(x) sech2(x).(3.13) Note that in our case V1(x) is an even function of xand V2(x) is an odd function. In our study we will assume ρ(x, t) = ˜ρ(y, t) where y(t) = x−q(t). That is, the functional form of ρwill be maintained if it is given a slight perturbation away from the origin. If it stays at the origin (q(t) = 0) and only changes its width and amplitude under perturbation, then we see that since ρ is an even function of yand V2(x) is an odd function of x, the mass is conserved. One can in a systematic fashion obtain the equations for the higher moments of hxnˆpmi, where ˆp=−i∂/∂x. It can be demonstrated that the four and six collective coordinate approximations we derive in this paper will satisfy a particular subset of four or six moment equations [35]. IV. COLLECTIVE COORDINATES The time dependent variational approximation relies on introducing a finite set of time-dependent real parameters in a trial wave function that one hopes captures the time evolution of a perturbed solution. By doing this one obtains a simplified set of ordinary differential equations for the collective coordinates in place of solving the full partial differential equation for the NLS equation. By judiciously choosing the collective coordinates, they can be simply related to the moments of xand ˆp=−i∂/∂x averaged over the density ρ(x, t). That is, we set ψ(x, t)7→ ˜ ψ[x, Q(t) ] (4.1) Q(t) = {Q1(t), Q2(t), . . . , Q2n(t)} ∈ R2n. The success of the method depends greatly on the choice of the the trial wave function ˜ ψ[x, Q(t) ]. The generalized Euler-Lagrange equations lead to Hamilton’s equations for the collective coordinates Q(t). Introducing the notation ∂µ≡∂/∂Qµ, the Lagrangian in terms of the collective coordinates is given by L(Q, ˙ Q) = πµ(Q)˙ Qµ−H(Q),(4.2) where πµ(Q) is defined by πµ(Q) = i 2Zdx{˜ ψ∗(x, Q) [ ∂µ˜ ψ(x, Q) ] (4.3) −[∂µ˜ ψ∗(x, Q) ] ˜ ψ(x, Q)}, and H(Q) is given by H(Q) =Zdxn|∂x˜ ψ(x, Q)|2−g 2|˜ ψ(x, Q)|4(4.4) +V1(x)|˜ ψ(x, Q)|2o. Similarly, in terms of the collective coordinates, the dissipation functional is given by F[Q, ˙ Q] = wµ(Q)˙ Qµ,(4.5) where wµ(Q) = iZdx V2(x){˜ ψ∗(x, Q) [ ∂µ˜ ψ(x, Q) ] (4.6) −[∂µ˜ ψ∗(x, Q) ] ˜ ψ(x, Q)}. The generalized Euler-Lagrange equations are ∂L ∂Qµ−d dt∂L ∂˙ Qµ=−∂F ∂˙ Qµ.(4.7) Setting vµ(Q) = ∂µH(Q), we find fµν(Q)˙ Qν=uµ(Q) = vµ(Q)−wµ(Q) (4.8) where fµν(Q) = ∂µπν(Q)−∂νπµ(Q) (4.9) 5 is an antisymmetric 2n×2nsymplectic matrix. If det{f(Q)} 6= 0, we can define an inverse as the contravariant matrix with upper indices, fµν(Q)fνσ(Q) = δµ σ,(4.10) in which case the equations of motion (4.8) can be put in the symplectic form: ˙ Qµ=fµν(Q)uν(Q).(4.11) Poisson brackets are defined using fµν(Q). If A(Q) and B(Q) are functions of Q, Poisson brackets are defined by {A(Q), B(Q)}= (∂µA(Q)) fµν (Q) (∂νB(Q)) .(4.12) In particular, {Qµ, Qν}=fµν(Q).(4.13) It is easy to show that fµν (x) satisfies Bianchi’s identity. This means that definition (4.12) satisfies Jacobi’s identity, as required for symplectic variables. The rate of energy loss is expressed as dH(Q) dt=−vµ(Q)fµν(Q)wν(Q), since fµν(Q) is an antisymmetric tensor. V. FOUR PARAMETER TRIAL WAVE FUNCTION Let us first look at the four parameter trial wave function that we have successfully used to study the effect of weak complex external potentials on the exact solution of the NLS equation in the absence of that potential. That is we will choose: ˜ ψ(x, t) = A0β(t) sech[ β(t)y(x, t) ] ei˜ φ(x,t),(5.1) where A0is the amplitude of the exact solution in the presence of the external potential (2.10) and is a funcition of m, b, g, and ˜ φ(x, t) = −θ(t) + p(t)y(x, t) + φ(x).(5.2) Here φ(x) is given by Eq. (2.11) and we have put y(x, t) = x−q(t). The four variational parameters are labeled by Qµ={q(t), p(t), β(t), θ(t)}.(5.3) The derivatives of ˜ ψ(x, t) with respect to tand xare given by ˜ ψt(x, t) = A0{˙ βsech(βy) (5.4a) −βsech(βy) tanh(βy) [ ˙ βy −˙qβ ] + i βsech(β y ) [ −˙ θ+ ˙p y −p˙q]}ei˜ φ(x,t), ˜ ψx(x, t) = A0β{ −βsech(βy) tanh(βy) (5.4b) + i sech(βy) [ p+ (2/3) m b sech(x) ] }ei˜ φ(x,t), where we have used (2.11b). Then the density and current is given by ρ(x, t) = A2 0β2sech2(βy),(5.5a) j(x, t)=2ρ(x, t) [ p+ (2/3) m b sech(x) ] .(5.5b) The time dependent mass, M(t) which is a normalization factor, is given by M(t) =Zdx ρ(x, t) = 2 A2 0β(t),(5.6) and the Lagrangian and dissipation function are given by, L=i 2Zdx[ψ∗ψt−ψ∗ tψ]−H[ψ, ψ∗],(5.7a) H=Zdx[|ψx|2−g|ψ|4/2 + V1(x)|ψ|2],(5.7b) F= iZdx V2(x) [ ψ∗ψt−ψ∗ tψ].(5.7c) The generalized Euler-Lagrange equations are δL δψ∗−∂t δL δψ∗ t =−δF δψ∗ t ,(5.8a) δL δψ −∂t δL δψt =−δF δψt .(5.8b) For the trial wave function of Eq. (5.1), we find L0[Q]≡i 2Zdx[˜ ψ∗˜ ψt−˜ ψ∗ t˜ ψ] (5.9) = 2A2 0β(˙ θ+p˙q)≡πµ(Q)˙ Qµ, where πq= 2A2 0β p, πp= 0, πβ= 0, πθ= 2A2 0β . (5.10) The only partial derivatives of πµ(Q) that survive are: ∂pπq= 2A2 0β, ∂βπq= 2A2 0p, ∂βπθ= 2A2 0.(5.11) So the symplectic matrix and its inverse are given by fµν(Q)=2A2 0   0−β−p0 β0 0 0 p0 0 1 0 0 −1 0   ,(5.12) fµν(Q) = 1 2A2 0β   0 1 0 0 −1 0 0 p 0 0 0 −β 0−p β 0   . From the Hamiltonian (5.7b) and our choice of trial wave function we find that H(Q) = A2 0β{(2/3) β2+ 2 p2+ (4/3) p m b β I1(β, q) −[b2+m2−(4/9) m2b2−1/4 ] β I2(β, q)}(5.13) −(2/3) g A4 0β3, 6 where I1(β, q) and I2(β, q) are given in Appendix A. Then defining vµ=∂µH(Q), we find vq=−A2 0β[ (4/3) p m b β f1(β, q) (5.14a) −2 [ b2+m2−(4/9) m2b2−1/4 ] β f6(β, q) ] , vp=A2 0β[ 4p+ (4/3) m b β I1(β, q) ] ,(5.14b) vβ=A2 0β{2β+ 2 p2/β (5.14c) + (8/3) p m b [I1(β, q)−β f10(β, q) ] −2 [ b2+m2−(4/9) m2b2−1/4 ] ×[I2(β, q)−β f7(β, q) ] −2g A2 0β} vθ= 0 ,(5.14d) where the fi(β, q) are given in Appendix A. From (5.7c), the dissipation function is given by F[Q, ˙ Q] = wµ(Q)˙ Qµ,(5.15) where wq=−4m b A2 0β2p f1(β, q),(5.16a) wp= 4 m b A2 0β2f2(β, q),(5.16b) wβ= 0 ,(5.16c) wθ=−4m b A2 0β2f1(β, q).(5.16d) Here f1(β, q) and f2(β, q) are given in Appendix A. In terms of the vector uµ(Q) = vµ(Q)−wµ(Q), Hamilton’s equations for the variational parameters are ˙ Qµ=fµν(Q)uν(Q),(5.17) which gives ˙q= 2 p+ (2/3) m b β I1(β, q)−2m b β f2(β, q),(5.18a) ˙p= (2/3) p m b β f1(β, q) (5.18b) −[b2+m2−(4/9) m2b2−1/4 ] β f6(β, q) ˙ β=−2β2m b f1(β, q).(5.18c) The equation for ˙ θis not needed for the evolution of the set of equations given in (5.18). For m= 0, the equations reduce to: ˙q= 2 p , (5.19a) ˙p=−[b2−1/4 ] β f6(β, q),(5.19b) ˙ β= 0 .(5.19c) So in this case, β= 1 and is fixed. This is because the normalization must be conserved. Equations (5.19) then reduce to: ¨q+ 2 [ b2−1/4 ] f6(1, q)=0.(5.20) A. Small Oscillation equations Using the expansions found in Appendix A we obtain for the small oscillation equations (we set q=δq,p=δp, 0.6 0.8 1.0 1.2 1.4 b 5 10 15 20 25 T FIG. 1. Period as a function of bfor m= 0 (upper curve) and m= 1 (lower curve) for 4 CC approximation. and β= 1 + δβ with δQµassumed small), δ˙q=π 72 9π2−64bm δβ + 2 δp , (5.21a) δ˙p=−8 15b2+m2−(4/9) b2m2−1/4δq , (5.21b) δ˙ β=−π 2mb δq . (5.21c) Thus we obtain for ¨q δ¨q+ω2(b, m)δq = 0 ,(5.22) where ω2(b, m) = π2 144 ( 9π2−64 ) b2m2(5.23) +16 15 b2+m2−(4/9) b2m2−1/4. The period T= 2π/ω(b, m) for m= 0 and m= 1 is shown in Fig. 1. VI. SIX PARAMETER ANSATZ One expects that when one increases the number of CC’s the accuracy of the variational approximation increases. For the six parameter Ansatz we will introduce a “chirp” term [39] Λ(t) which is conjugate to the width parameter β(t). That is we will assume: ˜ ψ(x, t) = A(t) sech[ β(t)y(x, t) ] ,ei˜ φ(x,t)(6.1) where ˜ φ(x, t) = −θ(t)+p(t)y(x, t)+Λ(t)y(x, t)2+φ(x).(6.2) Here φ(x) is given by Eq. (2.11) and we have put y(x, t) = x−q(t). We find ρ(x, t) = |˜ ψ(x, t)|2=A2(t) sech2(βy),(6.3) 7 so that the mass becomes M(t) =Zdx ρ(x, t) = 2A2(t) β(t).(6.4) It will be useful to employ M(t) as a collective coordinate rather than A(t). The six time-dependent collective coordinates then are: Qµ(t) = {M(t), θ(t), q(t), p(t), β(t),Λ(t)}.(6.5) The parameters β(t) and Λ(t) are related to the two point correlation functions G2=h(x−q(t))2iand P2= h[x−q(t)]ˆp+ ˆp[x−q(t)]iwhere h(·)i= ∞ Z −∞ (·)|ψ(x, t)|2dx. ∞ Z −∞ |ψ(x, t)|2dx . (6.6) Thus we find G2=π2/(12β2), and P2=i 2Zdx[x−q(t)] [ ψ∗ψx−ψ∗ xψ]/M(t) (6.7) =π2Λ 3β2+2 3bmI3(β, q) M(t), where I3is given in Appendix A. We see that P2is directly related to Λ when the potential is real. From the formalism given in Sec. IV, the equations of motion for the collective coordinates follow. For the kinetic term in the Lagrangian, we find πM= 0, πθ=M, πq=Mp, πp= 0 (6.8) πβ= 0, πΛ=−Mπ2 12β2, and the only non-zero derivatives are then ∂Mπθ= 1, ∂Mπq=p, ∂pπq=M(6.9) ∂MπΛ=−π2 12β2, ∂βπΛ=Mπ2 6β3. The antisymmetric symplectic tensor is then given by fµν(Q) =        0 1 p0 0 −π2/(12β2) −1 0 0 0 0 0 −p0 0 −M0 0 0 0 M0 0 0 0 0 0 0 0 Mπ2/(6β3) π2/(12β2) 0 0 0 −Mπ2/(6β3) 0        .(6.10) Since det{fij(Q)}=M4π4/(36β6) and is non-zero, the inverse is given by fµν(Q) =        0−1 0 0 0 0 1 0 0 −p/M β/(2M) 0 0 0 0 1/M 0 0 0p/M −1/M 0 0 0 0−β/(2M) 0 0 0 −6β3/(π2M) 0 0 0 0 6β3/(π2M) 0        .(6.11) For the dissipation functional, we obtain F(Q, ˙ Q)=2MmbβZdysech2(βy) sech(y+q) tanh(y+q) [ −˙ θ+ ˙py −p˙q+˙ Λy2−2yΛ ˙q],(6.12) which gives wM= 0, wθ=−2Mmbβ f1(β, q), wq=−2M mbβ [p f1(β, q) + 2Λf2(β, q) ],(6.13) wp= 2M mbβ f2(β, q), wβ= 0, wΛ= 2M mbβ f3(β, q). For H(Q), using the 6-parameter Ansatz we now obtain H(Q) M=p2+β2 3+π2Λ2 3β2+2β 3pbm I1(β, q) + 4β 3bmΛI3(β, q) (6.14) −gMβ 6−β 2hb2+m2−1 4−4 9b2m2iI2(β, q). 8 All the integrals are defined in Appendix A. For vµ(Q) = ∂µH(Q) we obtain vM=p2+β2 3+π2Λ2 3β2+2β 3pbm I1(β, q) + 4β 3bmΛI3(β, q) (6.15a) −gMβ 3−hb2+m2−1 4−4 9b2m2iβ 2I2(β, q), vθ= 0 ,(6.15b) vq=−2β 3M pbm f1(β, q)−4β 3M bmΛf2(β, q) (6.15c) +Mhb2+m2−1 4−4 9b2m2iβ f6(β, q), vp= 2Mp +2Mβ 3bm I1(β, q),(6.15d) vβ=2Mβ 3−2Mπ2Λ2 3β3+2M 3pbm I1(β, q) + 4M 3bmΛI3(β, q) (6.15e) −gM2 6−M 2hb2+m2−1 4−4 9b2m2iI2(β, q)−4Mβ 3pbm f10(β, q) −8Mβ 3bmΛf9(β, q) + hb2+m2−1 4−4 9b2m2iMβ f7(β, q), vΛ=2π2MΛ 3β2+4βM 3bm I3(β, q).(6.15f) The symplectic equations of motion are ˙ Qµ=fµν(Q)uν(Q),(6.16) from which we find: ˙ M=−2M mbβf1(β, q),(6.17a) ˙ θ=−p2+2 3β2−5 12 gβM +1 3mbpβ I1(β, q)+2mbβΛI3(β, q) (6.17b) + 2 mbpβ f2(β, q)−2 3mbpβ2f10(β, q)−4 3mbβ2Λf9(β, q) −1 4hb2+m2−1 4−4 9b2m2iβ[ 3 I2(β, q)−2β f7(β, q) ] ˙q= 2p+2β 3mb I1(β, q)−2mbβ f2(β, q),(6.17c) ˙p=2 3mbβ p f1(β, q)−8 3mbβ Λf2(β, q)−hb2+m2−1 4−4 9m2b2iβ f6(β, q) (6.17d) ˙ β=−mb β2f1(β, q)−4βΛ−8β4 π2mb I3(β, q) + 12β4 π2mb f3(β, q),(6.17e) ˙ Λ = −4Λ2+4β4 π2+4 π2β3pmb I1(β, q) + 8 π2β3Λmb I3(β, q) (6.17f) −gβ3M π2−6β3 π2hb2+m2−1 4−4 9b2m2if8(β, q) −8β4 π2bm p f10(β, q)−16β4 π2bmΛf9(β, q). In Eq. (6.17f), we use the identity (A10). Here M(t) is a dynamic variable. In order for the variational trial wave function to match the exact solution at t= 0, the initial conditions are: q0= 0, p0= 0, β0= 1,Λ0= 0, θ0=−t(6.18) gM0=(4b2−9)(4m2−9) 18 . 9 As a check, the right-hand-sides of Eqs. (6.17) vanish [except for ˙ θ(0) = −1] at these initial values, which guarantees that the exact solution is stationary. For non-zero values of q0and/or β0, the values of p0and Λ0are sometimes fixed by setting ˙q0= 0 and ˙ β0= 0, and solving Eqs. (6.17c), and (6.17e) for p0and Λ0, which gives: p0=1 2h˙q0−2 3mbβ0I1(β0, q0) (6.19a) + 2 mbβ0f2(β0, q0)i, Λ0=1 4β0h−˙ β0−mbβ2 0f1(β0, q0) (6.19b) −8 π2mbβ4I3(β0, q0) + 12 π2mbβ4f3(β0, q0)i. When m= 0, the external potential is real and ˙ M= 0. The stability of the solutions to this equation for arbitrary κand for repulsive and attractive potential V1as well as positive and negative gwas studied using a variety of methods, and the stability properties and small oscillation frequencies for q, p, β, Λ were determined in Ref. [38]. For that problem when we set κ= 1 and m= 0, our equations simplify to ˙q= 2 p , (6.20) ˙ β=−4βΛ, ˙ Λ = −4Λ2+4β4 π2−gβ3M π2−6β3 π2hb2−1 4if8(β, q) which agrees with the results in Ref. [38] once we use the fact that f3[G, q, γ] in that paper is just β2f8(β, q) here. At m= 0 the small oscillation equations for βand q decouple. Using the expansions of the integrals found in Appendix A, we find that the small oscillation equations are: δ˙q= 2δp (6.21) δ˙p=−8 15(b2−1/4) δq so that δ¨q+ω2 qδq = 0 ,(6.22) ω2 q=16 15(b2−1/4) . This agrees with the result from the 4-parameter Ansatz. However, we get a different frequency for the βoscillation, δ˙ β=−4δΛ,(6.23) δ˙ Λ = h4b2 15 +4 π2−1 15 iδβ , so that δ¨ β+ω2 βδβ = 0 ,(6.24) ω2 β= 4 h4b2 15 +4 π2−1 15 i. Plots of ω2 qand ω2 βfor m= 0 are shown in Fig. 2(a). VII. LINEAR RESPONSE RESULTS FOR THE SIX CC APPROXIMATION We linearize the set of equations given in (6.17) by expanding the equations about the exact solutions, Qµ= Qµ 0+δQµkeeping only the first order terms. Note that Qµ 0are given in Eqs. (6.18). Using the expansions of Appendix ??, we find δ˙ M=−π 2mb M0δq , (7.1a) δ˙ θ=−5 12 gδM +7π 18 mb δp (7.1b) +1 3h1 + 2π2 15 −1 2+π2 30 gM0iδβ , δ˙q=π 72 9π2−64 mb δβ + 2 δp , (7.1c) δ˙p=4 15 [gM0−4 ] δq −4π 9mb δΛ,(7.1d) δ˙ β=hπ 2−20 3πimb δq −4δΛ,(7.1e) δ˙ Λ = 2bm 3πδp −g π2δM (7.1f) +2 15h−gM0+30 π2+ 4 iδβ , where we have used the relation, b2+m2−4 9b2m2−1 4= 2 −1 2gM0.(7.2) Equations (7.1) are written as δ˙ Qµ=Mµν(Q0)δQν,(7.3) from which we find: δ¨ Qµ+Wµν(Q0)δQν= 0 (7.4) Wµν(Q0) = −Mµσ(Q0)Mσν(Q0). Here Wµν(Q0) is Hermitian. The square of the linearized oscillation frequencies ω2are given by the eigenvalues of Wµν(Q0). One can show that the matrix Wµν(Q0) can be split into two blocks, one of them coupling (δq, δΛ, δθ), the other coupling (δp, δβ, δM). Both of these blocks give identical eigenvalues, a zero eigenvalue and two non-zero eigenvalues. For example, using Eqs. (7.1), we find δ¨q−[A δq +B δΛ ] = 0 ,(7.5a) δ¨ Λ−[D δq +E δΛ ] = 0 ,(7.5b) 16 That is, I2(β, q)−2βf7(β, q) = 2 f8(β, q).(A10) We use this identity in the ˙ Λ equation, (6.17f). Next, we now consider the expansion of the integrals and find to first order: I1(1 + δβ, δq) = π 2−π 3δβ , (A11a) I2(1 + δβ, δq) = 4 3−2 3δβ , (A11b) I3(1 + δβ, δq) = −π 6δq , (A11c) f1(1 + δβ, δq) = π 4δq , (A11d) f2(1 + δβ, δq) = π 6+π 48(16 −3π2)δβ , (A11e) f3(1 + δβ, δq) = π 48(−32 + 3π2)δq , (A11f) f6(1 + δβ, δq) = 8 15 δq , (A11g) f7(1 + δβ, δq) = 1 3+2 45(−15 + π2)δβ , (A11h) f8(1 + δβ, δq) = 1 3−2π2 45 δβ , (A11i) f9(1 + δβ, δq) = π 96(16 −3π2)δq , (A11j) f10(1 + δβ, δq) = π 6+π 32(−16 + π2)δβ . (A11k) Appendix B: Generalized traveling wave method This method was named and used in a paper by Quintero, Mertens and Bishop [35]. We will show here that it is an alternative way to obtain Eq. (6.16) for the rate of change of the collective coordinates. The authors substitute the trial wave function directly into Schr¨odinger’s equation. This gives i˙ Qν∂ν˜ ψ(x, Q) + ˜ ψxx(x, Q) + g|˜ ψ(x, Q)|2κ˜ ψ(x, Q) = [ V1(x)+iV2(x) ] ˜ ψ(x, Q),(B1a) −i˙ Qν∂ν˜ ψ∗(x, Q) + ˜ ψxx(x, Q) + g|˜ ψ(x, Q)|2κ˜ ψ∗(x, Q) = [ V1(x)−iV2(x) ] ˜ ψ∗(x, Q).(B1b) Multiply (B1a) by ∂µ˜ ψ∗(x, Q) and (B1b) by ∂µ˜ ψ(x, Q) and add them to give i{[∂µ˜ ψ∗] [∂ν˜ ψ]−[∂ν˜ ψ∗] [∂µ˜ ψ]}˙ Qν+ [∂µ˜ ψ∗]˜ ψxx (B2) + [∂µ˜ ψ]˜ ψ∗ xx +{g|˜ ψ|2κ−V1(x)} { [∂µ˜ ψ∗]˜ ψ+ [∂µ˜ ψ]˜ ψ∗} = i V2(x){[∂µ˜ ψ∗]˜ ψ−[∂µ˜ ψ]˜ ψ∗}. Integrating (B2) over xand the second term by parts gives Iµν(Q)˙ Qν=∂µH(Q) + Rµ(Q),(B3) where Iµν(Q) = i Zdx[∂µ˜ ψ∗] [ ∂ν˜ ψ]−[∂ν˜ ψ∗] [ ∂µ˜ ψ],(B4a) H(Q) =Zdx|∂x˜ ψ|2−g|˜ ψ|2κ+2 κ+ 1 +V1(x)|˜ ψ|2,(B4b) Rµ(Q)=iZdx V2(x)[∂µ˜ ψ∗]˜ ψ−˜ ψ∗[∂µ˜ ψ].(B4c) Here we have interchanged µ↔νin the definition of Iµν(Q) from their Eq. (6) [35]. So we see that Rµ(Q)≡ −wµ(Q) and we find that 17 fµν(Q) = ∂µπν(Q)−∂νπµ(Q) (B5) =i 2Zdx{[∂µ˜ ψ∗] [ ∂ν˜ ψ] + ˜ ψ∗[∂µ∂ν˜ ψ]−[∂µ∂ν˜ ψ∗]˜ ψ−[∂ν˜ ψ∗] [ ∂µ˜ ψ]−[∂ν˜ ψ∗] [ ∂µ˜ ψ]−˜ ψ∗[∂ν∂µ˜ ψ] + [ ∂ν∂µ˜ ψ∗]˜ ψ+ [ ∂µ˜ ψ∗] [ ∂ν˜ ψ]}, = iZdx{[∂µ˜ ψ∗] [ ∂ν˜ ψ]−[∂ν˜ ψ∗] [ ∂µ˜ ψ]}=Iµν(Q). In the notation used in the variational method, Eq. (B3) becomes fµν(Q)˙ Qν=uµ(Q)−wµ(Q) = vµ(Q).(B6) So the generalized traveling wave approximation is identical to the variational method. The authors of Ref. [35] proved this in another way in Sec. III of their paper for a simpler dissipative system. [1] C. M. Bender, Rep. Prog. Phys. 70, 947 (2007). [2] H. Geyer, D. Heiss, and M. Znojil, J. Phys. A: Math. Gen. 39 (2006). [3] A. Fring, H. Jones, and M. Znojil, J. Phys. A: Math. Theor. 41, 240301 (2008). [4] C. M. Bender, A. Fring, U. G¨unther, and H. Jones, J. Phys. A: Math. Theor. 45, 440301 (2012). [5] K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani, Int. J. Theor. Phys. 50, 1019 (2011). [6] A. Ruschhaupt, F. Delgado, and J. G. Muga, J. Phys. A: Math. Gen. 38, L171 (2005). [7] K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani, Phys. Rev. Lett. 100, 103904 (2008). [8] S. Klaiman, U. G¨unther, and N. Moiseyev, Phys. Rev. Lett. 101, 080402 (2008). [9] S. Longhi, Phys. Rev. Lett. 103, 123601 (2009). [10] S. Longhi, Phys. Rev. B 80, 235102 (2009). [11] S. Longhi, Phys. Rev. A 81, 022102 (2010). [12] C. E. Ruter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, and D. Kip, Nat. Phys. 6, 192 (2010). [13] A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G. A. Siviloglou, and D. N. Christodoulides, Phys. Rev. Lett. 103, 093902 (2009). [14] A. Regensburger, C. Bersch, M.-A. Miri, G. Onishchukov, D. N. Christodoulides, and U. Peschel, Nature 488, 167 (2012). [15] J. Schindler, A. Li, M. C. Zheng, F. M. Ellis, and T. Kottos, Phys. Rev. A 84, 040101 (2011). [16] J. Schindler, Z. Lin, J. M. Lee, H. Ramezani, F. M. Ellis, and T. Kottos, J. Phys. A: Math. Theor. 45, 444029 (2012). [17] C. M. Bender, B. K. Berntson, D. Parker, and E. Samuel, Am. J. Phys. 81, 173 (2013). [18] B. Peng, S. K. Ozdemir, F. Lei, F. Monifi, M. Gianfreda, G. L. Long, S. Fan, F. Nori, C. M. Bender, and L. Yang, Nat. Phys. 10, 394 (2014). [19] M.-A. Miri, M. Heinrich, R. El-Ganainy, and D. N. Christodoulides, Phys. Rev. Lett. 110, 233902 (2013). [20] M. Heinrich, M.-A. Miri, S. St¨utzer, R. El-Ganainy, S. Nolte, A. Szameit, and D. N. Christodoulides, Nature Communications 5, 3698 EP (2014). [21] L. ´ E. Gendenshte˘ın and I. V. Krive, Soviet Physics Uspekhi 28, 645 (1985). [22] F. Cooper, A. Khare, and U. Sukhatme, Phys. Rept. 251, 267 (1995). [23] B. Bagchi and R. Roychoudhury, J. Phys. A: Math. Gen. 33, L1 (2000). [24] B. Bagchi, S. Mallik, and C. Quesne, Int. J. Mod. Phys. A16, 2859 (2001). [25] B. Bagchi and C. Quesne, Phys. Lett. A 273, 285 (2000). [26] Z. Ahmed, Phys. Lett. A 282, 343 (2001). [27] B. Midya, Phys. Rev. A 89, 032116 (2014). [28] P. G. Kevrekidis, J. Cuevas-Maraver, A. Saxena, F. Cooper, and A. Khare, Phys. Rev. E 92, 042901 (2015). [29] F. Cooper, A. Khare, and U. Sukhatme, Supersymmetry in Quantum Mechanics (World Scientific, Singapore; River Edge, N.J., 2001). [30] F. Cooper, A. Khare, A. Comech, B. Mihaila, J. F. Dawson, and A. Saxena, J. Phys. A: Math. Theor. 50, 015391 (2017). [31] A. Comech, arXiv:1203.3859 (2012). [32] N. G. Vakhitov and A. A. Kolokolov, Radiophys. Quantum Electron 16, 783 (1973). [33] F. G. Mertens, F. Cooper, E. Ar´evalo, A. Khare, A. Saxena, and A. R. Bishop, Phys. Rev. E 94, 032213 (2016). [34] F. Cooper, A. Khare, N. R. Quintero, F. G. Mertens, and A. Saxena, Phys. Rev. E 85, 046607 (2012). [35] N. R. Quintero, F. G. Mertens, and A. R. Bishop, Phys. Rev. E 82, 016606 (2010). [36] G. Poschl and E. Teller, Z. Phys. 83, 143 (1933). [37] L. D. Landau and E. M. Lifshitz, Quantum Mechnics (Nauka, Moscow, 1989). [38] J. F. Dawson, F. Cooper, A. Khare, B. Mihaila, E. Arevalo, R. Lan, A. Comech, and A. Saxena, “Stability of new exact solutions of the nonlinear Schr¨odinger equation in a Poschl-Teller external potential,” arXiv:1705.07253v1 (2017). [39] F. Cooper, S.-Y. Pi, and P. N. Stancioff, Phys. Rev. D 18 34, 3831 (1986). [40] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in FORTRAN: The Art of Scientific Computing (Cambridge University Press, Cambridge, 1992).