Response of exact solutions of the nonlinear Schrodinger equation to small perturbations in a class of complex external potentials having supersymmetry and parity-time symmetry
Fred Cooper, John F. Dawson, Franz G. Mertens, Edward Arevalo, Niurka R. Quintero, Bogdan Mihaila, Avinash Khare, Avadh Saxena
LLA-UR-17-25332
Response of exact solutions of the nonlinear Schr¨odinger equation to smallperturbations in a class of complex external potentials having supersymmetry andparity-time symmetry
Fred Cooper,
1, 2, ∗ John F. Dawson, † Franz G. Mertens, ‡ Edward Ar´evalo, § NiurkaR. Quintero,
6, 7, ¶ Bogdan Mihaila, ∗∗ Avinash Khare, †† and Avadh Saxena ‡‡ The Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA Department of Physics, University of New Hampshire, Durham, NH 03824, USA Physikalisches Institut, Universit¨at Bayreuth, D-95440 Bayreuth, Germany Pontifical Catholic University of Chile, Instituto de F´ısica, Santiago, Regi´on Metropolitana, Chile Instituto Carlos I de Fisica Te´orica y Computacional,Universidad de Granada, E-18015, Granada, Spain Departamento de Fisica Aplicada I, E.P.S. Universidad de Sevilla, 41011 Sevilla, Spain Physics Division, National Science Foundation, Arlington, VA 22230, USA Physics 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¨odingerequation in 1+1 dimensions in an external complex potential derivable from a parity-time symmetricsuperpotential that was considered earlier [Phys. Rev. E 92, 042901 (2015)]. In particular weconsider the nonlinear partial differential equation { i ∂ t + ∂ x + g | ψ ( x, t ) | − V + ( x ) } ψ ( x, t ) = 0, where V + ( x ) = (cid:0) − b − m + 1 / (cid:1) sech ( x ) − i m b sech( x ) tanh( x ) represents the complex potential. Herewe study the perturbations as a function of b and m using a variational approximation based ona dissipation functional formalism. We compare the result of this variational approach with directnumerical simulation of the equations. We find that the variational approximation works quite wellat small and moderate values of the parameter bm which controls the strength of the imaginarypart of the potential. We also show that the dissipation functional formalism is equivalent to thegeneralized 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 applica-tions on the one hand, as well as its mathematical struc-ture on the other, have drawn considerable attention fromboth the physics and the mathematics community. Theoriginal proposal of Bender and his collaborators [1–4]towards the study of such systems was made as an al-ternative to the postulate of Hermiticity in quantum me-chanics. Keeping in perspective the formal similarity ofthe Schr¨odinger equation with Maxwell’s equations in theparaxial approximation, it was realized that such PT in-variant systems can in fact be experimentally realized inoptics [5–14]. Subsequently, these efforts motivated ex-periments in several other areas including PT invariantelectronic circuits [15, 16], mechanical circuits [17], andwhispering-gallery microcavities [18].Concurrently, the notion of supersymmetry (SUSY)originally espoused in high-energy physics has also been ∗ [email protected] † [email protected] ‡ [email protected] § earevalo@fis.puc.cl ¶ [email protected] ∗∗ [email protected] †† [email protected] ‡‡ [email protected] realized in optics [19, 20]. The key idea is that from agiven potential one can obtain a SUSY partner potentialwith both potentials possessing the same spectrum, ex-cept possibly for one eigenvalue [21, 22]. Therefore, aninterplay of SUSY with PT symmetry is expected to bequite rich and is indeed very useful in achieving trans-parent as well as one-way reflectionless complex opticalpotentials [23–27].A previous paper [28] explored the interplay between PT symmetry, SUSY and nonlinearity. That paper de-rived 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, theyconsidered the nonlinear partial differential equation (cid:8) i ∂ t + ∂ x − V ± ( x ) + g | ψ ( x, t ) | κ (cid:9) ψ ( x, t ) = 0 , (1.1)for arbitrary nonlinearity parameter κ , with V ± ( x ) = W ( x ) ∓ W (cid:48) ( x ) − ( m − / , (1.2)and the partner potentials arise from the superpotential W ( x ) = ( m − /
2) tanh x − ib sech x , (1.3)giving rise to V + ( x ) = (cid:0) − b − m + 1 / (cid:1) sech ( x ) (1.4a) − i m b sech( x ) tanh( x ) ,V − ( x ) = (cid:0) − b − ( m − + 1 / (cid:1) sech ( x ) (1.4b) − i ( m − b sech( x ) tanh( x ) . a r X i v : . [ n li n . PS ] J u l For m = 1, the complex potential V + ( x ) has the samespectrum, apart from the ground state, as the real po-tential V − ( x ) and this fact was used in the numericalstudy of the stability of the bound state solutions of theNLS equation in the presence of V + ( x ) (see Ref. [28]).In a recent complementary study [30] of this system ofnonlinear Schr¨odinger equations in PT symmetric SUSYexternal potentials, the stability properties of the boundstate solutions of NLS equation in the presence of theexternal real SUSY partner potential V − ( x ) were inves-tigated. The stability regime of these solutions, whichdepended on the parameters ( b, κ ), was compared to thestability regime of the related solitary wave solutions tothe NLS equation in the absence of the external poten-tial. Because the NLS equation in the presence of V − ( x )is a Hamiltonian dynamical system, in Ref. [30] they wereable to use several variational methods to study the sta-bility of the solutions when they undergo certain smalldeformations, and showed that these variational meth-ods agreed with a linear stability analysis based on theVakhitov-Kolokolov (V-K) stability criterion [31, 32] aswell as numerical simulations that have recently been per-formed.In Ref. [28] we determined the exact solutions of theequation for m = 1 for V + ( x ), which was complex. Westudied numerically the stability properties of these so-lutions using linear stability analysis. We found someunusual results for the stability which depended on thevalue of b . What was found for m = 1 (and κ = 1) wasthat the eigenvalues of the linear stability matrix becamecomplex for 0 . < b < . bm . In this paper wefocus on the external potential V + ( x ) which is symmetricin b ↔ m .Here we will compare the numerical simulations withthe results of our collective coordinate (CC) approxima-tion. We will also look at the linear stability analysisthat arises from studying the linearization of the CC or-dinary differential equations (ODEs). For the case of areal external potential, studying the eigenvalues of thisreduced stability analysis predicted the correct stabilityregime [34].This paper is structured as follows. In Sec. II we re-view the non-hermitian SUSY model that we studied in Ref. [28] and add the self-interactions of the NLS equa-tion to the linear model. In Sec. III we give some of theexact low order moment equations for this problem. InSec. IV we introduce our collective coordinate approach,whereas in Sec. V we use a four parameter trial wavefunction that we considered in an earlier study of solitonbehavior in complex periodic external potentials, and de-rive equations for the four CC’s. In Sec. VI we expandthe number of CC’s to six and derive equations for thesix CC’s. In Sec. VII we study the linear response theoryof the six CC approximation. In Sec. VIII we present ournumerical strategy for solving the NLS equation startingfrom a perturbed exact solution. In Sec. IX we comparethe four and six CC approximations with direct numer-ical simulations. In Sec. X we present our main conclu-sions. Finally in Appendix A we provide the definitionsof various integrals and in Appendix B we show that forthis problem our variational approach is equivalent to thegeneralized traveling wave method [35]. II. NLS EQUATION IN THE PRESENCE OF ANON-HERMITIAN SUPERSYMMETRICEXTERNAL POTENTIAL
We were interested in studying the NLS equation inthe presence of a complex external potential and wereintrigued by the fact that as a result of PT symmetry,there existed complex potentials whose SUSY partnerswere real and had explicitly known spectra of boundstates. This led us to study the external potential de-fined by the PT symmetric SUSY superpotential W ( x )given by Eq. (1.3). This superpotential gives rise to su-persymmetric partner potentials given by Eqs. (1.4). Forthe case m = 1, V − ( x ) is the well known P¨oschl-Tellerpotential [36, 37]. The relevant bound state eigenvaluesassume an extremely simple form as E ( − ) n = −
14 [ 2 b − n − . (2.1)Such bound state eigenvalues only exist when n < b − /
2. We notice that for the ground state (n=0) to existrequires b > /
2. The existence of a first excited state(n=1) requires b > /
2. Here we consider the general V + ( x ) arising from the superpotential W ( x ) dependingon m, b as an external potential modifying the nonlinearSchr¨odinger equation. Rewriting the external potentialgiven in Eq. (1.4a) as, V + ( x ) = V ( x ) + i V ( x ) , (2.2)we have V ( x ) = − ( b + m − /
4) sech ( x ) , (2.3a) V ( x ) = − mb 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 ofthe exact solutions of the NLS equation in this externalpotential: { i ∂ t + ∂ x + g | ψ ( x, t ) | κ − [ V ( x ) + iV ( 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Γ = (cid:90) d t (cid:110) i2 (cid:90) d x [ ψ ∗ ψ t − ψψ ∗ t ] − H (cid:111) , (2.6a) H = (cid:90) d x (cid:110) | ψ x | − g | ψ | κ +2 κ + 1 + V ( x ) | ψ | (cid:111) , (2.6b) F = (cid:90) d t F = i (cid:90) d x d t V ( x ) [ ψ t ψ ∗ − ψ ∗ t ψ ] . (2.6c)Localized solutions to Eq. (2.4) exist for arbitrary val-ues of κ, m, b . Here we use ψ ( x, t ) to denote the exactsolution to the NLS equation in the external potential, ψ ( x, t ) = A sech /κ ( x ) e i[ Et + φ ( x ) ] (2.7)where φ ( x ) = 4 bmκκ + 2 tan − [tanh( x/
2) ] , (2.8)with E = 1 /κ , and gA κ = [4 b κ − ( κ + 2) ] [4 m κ − ( κ + 2) ]4 κ ( κ + 2) . (2.9)We notice when mb = 0 the potential is real and thatsolutions exist for m + b − / < ( κ + 1) /κ . There aretwo regimes where A is positive and so a solution existswhen m (cid:54) = 0. This form of the solution reflects the factthat the potential V + is invariant under the interchange m ↔ b .In a previous paper [38] we studied the stability ofthese solutions for m = 0 (real external potential) andfor arbitrary κ . In that paper, we also considered twoother cases where exact solutions exist. For the case of g = − V = 0, all thesolutions that were allowed were stable. Solutions alsoexist for V (cid:54) = 0 and are given by Eqs. (2.7) and (2.9)with g = −
1. For g = 1 and a repulsive real potentialwe found the solutions for V = 0 were translationallyunstable. Solutions again exist when V (cid:54) = 0 for thiscase. We will not discuss these solutions further here.Here we will confine ourselves to κ = 1 , g = 1 and anattractive external potential V and study the domain ofapplicability of the variational methods we have devel-oped previously to the case of increasing the dissipationby allowing m to vary. In particular for the case we willconcentrate on here ( κ = 1) we have that gA = (4 b − m − / , (2.10) so that when m < / b < / b + m > /
4. Note that gA is inde-pendent of g . For κ = 1 we have φ ( x ) = (4 m b/
3) tan − [ tanh( x/
2) ] , (2.11a) ∂ x φ ( x ) = (2 / m b sech( x ) . (2.11b) III. SOME GENERAL PROPERTIES OF THENLS EQUATION IN COMPLEX POTENTIALS
We are interested in solitary wave solutions that ap-proach zero exponentially at ±∞ . For these solutions wedefine the mass density ρ ( x, t ) = | ψ ( x, t ) | , and the massor norm M ( t ) as M ( t ) = (cid:90) d x ρ ( x, t ) = (cid:90) d x | ψ ( x, t ) | . (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 sub-tracting the complex conjugate of the resulting equation,we obtain ∂ρ ( x, t ) ∂t + ∂j ( x, t ) ∂x = 2 V ( x ) ρ ( x, t ) . (3.3)Integrating over space, and assuming that j (+ ∞ , t ) − j ( −∞ , t ) = 0, we findd M ( t )d t = 2 (cid:90) d x V ( x ) ρ ( x, t ) . (3.4)Note that M is conserved when V ( x ) = 0. If we insteadmultiply the NLS equation by ψ ∗ and add the complexconjugate of the resulting equation, we geti ( ψ ∗ ψ t − ψψ ∗ t ) (3.5)= − gρ − ψ ∗ ψ xx − ψψ ∗ xx + 2 V ( x ) ρ , which when we integrate over space, leads to the virialtheorem:i2 (cid:90) d x ( ψ ∗ ψ t − ψ ∗ t ψ ) − (cid:90) d x (cid:104) | ψ x | − g | ψ | (cid:105) (3.6)= (cid:90) d x V ( x ) | ψ | . The average position q ( t ) can be defined through the firstmoment of x as follows: M ( t ) = (cid:90) d x x ρ ( x, t ) = q ( t ) M ( t ) . (3.7)Multiplying the continuity equation (3.3) by x and inte-grating over all space we find:d M d t = 2 P ( t ) + 2 (cid:90) d x x V ( x ) ρ ( x, t ) , (3.8)where the momentum P ( t ) = 12 (cid:90) d x j ( x, t ) (3.9)= i2 (cid:90) d x [ ψ ∗ ( x, t ) ψ x ( x, t ) − ψ ∗ x ( x, t ) ψ ( x, t ) ] . Here, we assumed thatlim 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 finddd t [ M ( t ) q ( t ) ] = 2 P ( t ) + 2 (cid:90) d y y V ( y + q ( t )) ρ ( y, t )+ 2 q ( t ) (cid:90) d x V ( x ) ρ ( x − q ( t ) , t ) . We recognize the last term as q ( t ) d M ( t ) / d t , so that wefinally have: M ( t ) d q ( t )d t = 2 P ( t )+2 (cid:90) d y y V ( y + q ( t )) ρ ( y, t ) . (3.11)Taking the time derivative of the momentum P ( t ), usingthe equations of motion for ψ and ψ ∗ , and integrating byparts, we findd P ( t )d t = − (cid:90) d x ρ ( x, t ) ∂V ( x ) ∂x + (cid:90) d x j ( x, t ) V ( x ) . (3.12)Here ∂V ( x ) ∂x = 2 (cid:0) b + m − / (cid:1) tanh( x ) sech ( x ) . (3.13)Note that in our case V ( x ) is an even function of x and V ( x ) is an odd function. In our study we will assume ρ ( x, t ) = ˜ ρ ( y, t ) where y ( t ) = x − q ( t ). That is, thefunctional form of ρ will be maintained if it is given aslight perturbation away from the origin. If it stays atthe origin ( q ( t ) = 0) and only changes its width andamplitude under perturbation, then we see that since ρ is an even function of y and V ( x ) is an odd function of x ,the mass is conserved. One can in a systematic fashionobtain the equations for the higher moments of (cid:104) x n ˆ p m (cid:105) ,where ˆ p = − i∂/∂x . It can be demonstrated that the fourand six collective coordinate approximations we derive inthis paper will satisfy a particular subset of four or sixmoment equations [35]. IV. COLLECTIVE COORDINATES
The time dependent variational approximation relieson introducing a finite set of time-dependent real param-eters in a trial wave function that one hopes captures thetime evolution of a perturbed solution. By doing this oneobtains a simplified set of ordinary differential equationsfor the collective coordinates in place of solving the fullpartial differential equation for the NLS equation. Byjudiciously choosing the collective coordinates, they canbe simply related to the moments of x and ˆ p = − i∂/∂x averaged over the density ρ ( x, t ).That is, we set ψ ( x, t ) (cid:55)→ ˜ ψ [ x, Q ( t ) ] (4.1) Q ( t ) = { Q ( t ) , Q ( t ) , . . . , Q n ( t ) } ∈ R n . The success of the method depends greatly on the choiceof the the trial wave function ˜ ψ [ x, Q ( t ) ]. The general-ized Euler-Lagrange equations lead to Hamilton’s equa-tions for the collective coordinates Q ( t ). Introducing thenotation ∂ µ ≡ ∂/∂Q µ , the Lagrangian in terms of thecollective coordinates is given by L ( Q, ˙ Q ) = π µ ( Q ) ˙ Q µ − H ( Q ) , (4.2)where π µ ( Q ) is defined by π µ ( Q ) = i2 (cid:90) d x { ˜ ψ ∗ ( x, Q ) [ ∂ µ ˜ ψ ( x, Q ) ] (4.3) − [ ∂ µ ˜ ψ ∗ ( x, Q ) ] ˜ ψ ( x, Q ) } , and H ( Q ) is given by H ( Q ) = (cid:90) d x (cid:110) | ∂ x ˜ ψ ( x, Q ) | − g | ˜ ψ ( x, Q ) | (4.4)+ V ( x ) | ˜ ψ ( x, Q ) | (cid:111) . Similarly, in terms of the collective coordinates, the dis-sipation functional is given by F [ Q, ˙ Q ] = w µ ( Q ) ˙ Q µ , (4.5)where w µ ( Q ) = i (cid:90) d x V ( x ) { ˜ ψ ∗ ( x, Q ) [ ∂ µ ˜ ψ ( x, Q ) ] (4.6) − [ ∂ µ ˜ ψ ∗ ( x, Q ) ] ˜ ψ ( x, Q ) } . The generalized Euler-Lagrange equations are ∂L∂Q µ − dd t (cid:16) ∂L∂ ˙ Q µ (cid:17) = − ∂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)is an antisymmetric 2 n × n symplectic matrix. Ifdet { f ( Q ) } (cid:54) = 0, we can define an inverse as the contra-variant matrix with upper indices, f µν ( Q ) f νσ ( Q ) = δ µσ , (4.10)in which case the equations of motion (4.8) can be putin 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 iden-tity. This means that definition (4.12) satisfies Jacobi’sidentity, as required for symplectic variables. The rate ofenergy loss is expressed asd H ( Q )d t = − v µ ( Q ) f µν ( Q ) w ν ( Q ) , since f µν ( Q ) is an antisymmetric tensor. V. FOUR PARAMETER TRIAL WAVEFUNCTION
Let us first look at the four parameter trial wave func-tion that we have successfully used to study the effect ofweak complex external potentials on the exact solution ofthe NLS equation in the absence of that potential. Thatis we will choose:˜ ψ ( x, t ) = A β ( t ) sech[ β ( t ) y ( x, t ) ] e i ˜ φ ( x,t ) , (5.1)where A is the amplitude of the exact solution in thepresence of the external potential (2.10) and is a funcitionof 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 t and x aregiven by˜ ψ t ( x, t ) = A { ˙ β sech( βy ) (5.4a) − β sech( βy ) tanh( βy ) [ ˙ βy − ˙ qβ ]+ i β sech( β y ) [ − ˙ θ + ˙ p y − p ˙ q ] } e i ˜ φ ( x,t ) , ˜ ψ x ( x, t ) = A β { − β sech( βy ) tanh( βy ) (5.4b)+ i sech( βy ) [ p + (2 / m b sech( x ) ] } e i ˜ φ ( x,t ) , where we have used (2.11b). Then the density and cur-rent is given by ρ ( x, t ) = A β sech ( βy ) , (5.5a) j ( x, t ) = 2 ρ ( x, t ) [ p + (2 / m b sech( x ) ] . (5.5b)The time dependent mass, M ( t ) which is a normalizationfactor, is given by M ( t ) = (cid:90) d x ρ ( x, t ) = 2 A β ( t ) , (5.6)and the Lagrangian and dissipation function are givenby, L = i2 (cid:90) d x [ ψ ∗ ψ t − ψ ∗ t ψ ] − H [ ψ, ψ ∗ ] , (5.7a) H = (cid:90) d x [ | ψ x | − g | ψ | / V ( x ) | ψ | ] , (5.7b) F = i (cid:90) d x V ( 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 L [ Q ] ≡ i2 (cid:90) d x [ ˜ ψ ∗ ˜ ψ t − ˜ ψ ∗ t ˜ ψ ] (5.9)= 2 A β ( ˙ θ + p ˙ q ) ≡ π µ ( Q ) ˙ Q µ , where π q = 2 A β p, π p = 0 , π β = 0 , π θ = 2 A β . (5.10)The only partial derivatives of π µ ( Q ) that survive are: ∂ p π q = 2 A β, ∂ β π q = 2 A p, ∂ β π θ = 2 A . (5.11)So the symplectic matrix and its inverse are given by f µν ( Q ) = 2 A − β − p β p − , (5.12) f µν ( Q ) = 12 A β − p − β − p β . From the Hamiltonian (5.7b) and our choice of trial wavefunction we find that H ( Q ) = A β { (2 / β + 2 p + (4 / p m b β I ( β, q ) − [ b + m − (4 / m b − / β I ( β, q ) } (5.13) − (2 / g A β , where I ( β, q ) and I ( β, q ) are given in Appendix A. Thendefining v µ = ∂ µ H ( Q ), we find v q = − A β [ (4 / p m b β f ( β, q ) (5.14a) − b + m − (4 / m b − / β f ( β, q ) ] ,v p = A β [ 4 p + (4 / m b β I ( β, q ) ] , (5.14b) v β = A β { β + 2 p /β (5.14c)+ (8 / p m b [ I ( β, q ) − β f ( β, q ) ] − b + m − (4 / m b − / × [ I ( β, q ) − β f ( β, q ) ] − g A β } v θ = 0 , (5.14d)where the f i ( β, q ) are given in Appendix A. From (5.7c),the dissipation function is given by F [ Q, ˙ Q ] = w µ ( Q ) ˙ Q µ , (5.15)where w q = − m b A β p f ( β, q ) , (5.16a) w p = 4 m b A β f ( β, q ) , (5.16b) w β = 0 , (5.16c) w θ = − m b A β f ( β, q ) . (5.16d)Here f ( β, q ) and f ( β, q ) are given in Appendix A. Interms of the vector u µ ( Q ) = v µ ( Q ) − w µ ( Q ), Hamilton’sequations for the variational parameters are˙ Q µ = f µν ( Q ) u ν ( Q ) , (5.17)which gives˙ q = 2 p + (2 / m b β I ( β, q ) − m b β f ( β, q ) , (5.18a)˙ p = (2 / p m b β f ( β, q ) (5.18b) − [ b + m − (4 / m b − / β f ( β, q )˙ β = − β m b f ( β, q ) . (5.18c)The equation for ˙ θ is not needed for the evolution of theset of equations given in (5.18). For m = 0, the equationsreduce to: ˙ q = 2 p , (5.19a)˙ p = − [ b − / β f ( β, q ) , (5.19b)˙ β = 0 . (5.19c)So in this case, β = 1 and is fixed. This is because thenormalization must be conserved. Equations (5.19) thenreduce to: ¨ q + 2 [ b − / f (1 , q ) = 0 . (5.20) A. Small Oscillation equations
Using the expansions found in Appendix A we obtainfor the small oscillation equations (we set q = δq , p = δp , FIG. 1. Period as a function of b for m = 0 (upper curve)and m = 1 (lower curve) for 4 CC approximation. and β = 1 + δβ with δQ µ assumed small), δ ˙ q = π (cid:0) π − (cid:1) bm δβ + 2 δp , (5.21a) δ ˙ p = − (cid:0) b + m − (4 / b m − / (cid:1) δq , (5.21b) δ ˙ β = − π mb δq . (5.21c)Thus we obtain for ¨ qδ ¨ q + ω ( b, m ) δq = 0 , (5.22)where ω ( b, m ) = π
144 ( 9 π −
64 ) b m (5.23)+ 1615 (cid:0) b + m − (4 / b m − / (cid:1) . The period T = 2 π/ω ( b, m ) for m = 0 and m = 1 isshown in Fig. 1. VI. SIX PARAMETER ANSATZ
One expects that when one increases the number ofCC’s the accuracy of the variational approximation in-creases. For the six parameter Ansatz we will introducea “chirp” term [39] Λ( t ) which is conjugate to the widthparameter β ( t ). That is we will assume:˜ ψ ( x, t ) = A ( t ) sech[ β ( t ) y ( x, t ) ] , e i ˜ φ ( x,t ) (6.1)where˜ φ ( x, t ) = − θ ( t ) + p ( t ) y ( x, t ) + Λ( t ) y ( x, t ) + φ ( 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 ) | = A ( t ) sech ( βy ) , (6.3)so that the mass becomes M ( t ) = (cid:90) d x ρ ( x, t ) = 2 A ( t ) β ( t ) . (6.4)It will be useful to employ M ( t ) as a collective coordi-nate rather than A ( t ). The six time-dependent collectivecoordinates then are: Q µ ( t ) = { M ( t ) , θ ( t ) , q ( t ) , p ( t ) , β ( t ) , Λ( t ) } . (6.5)The parameters β ( t ) and Λ( t ) are related to the twopoint correlation functions G = (cid:104) ( x − q ( t )) (cid:105) and P = (cid:104) [ x − q ( t )]ˆ p + ˆ p [ x − q ( t )] (cid:105) where (cid:104) ( · ) (cid:105) = ∞ (cid:90) −∞ ( · ) | ψ ( x, t ) | d x (cid:46) ∞ (cid:90) −∞ | ψ ( x, t ) | d x . (6.6) Thus we find G = π / (12 β ), and P = i2 (cid:90) d x [ x − q ( t )] [ ψ ∗ ψ x − ψ ∗ x ψ ] /M ( t ) (6.7)= π Λ3 β + 23 bm I ( β, q ) M ( t ) , where I is given in Appendix A. We see that P is di-rectly related to Λ when the potential is real.From the formalism given in Sec. IV, the equationsof motion for the collective coordinates follow. For thekinetic term in the Lagrangian, we find π M = 0 , π θ = M, π q = M p, π p = 0 (6.8) π β = 0 , π Λ = − M π β , and the only non-zero derivatives are then ∂ M π θ = 1 , ∂ M π q = p, ∂ p π q = M (6.9) ∂ M π Λ = − π β , ∂ β π Λ = M π β . The antisymmetric symplectic tensor is then given by f µν ( Q ) = p − π / (12 β ) − − p − M M M π / (6 β ) π / (12 β ) 0 0 0 − M π / (6 β ) 0 . (6.10)Since det { f ij ( Q ) } = M π / (36 β ) and is non-zero, the inverse is given by f µν ( Q ) = − − p/M β/ (2 M ) 00 0 0 1 /M p/M − /M − β/ (2 M ) 0 0 0 − β / ( π M )0 0 0 0 6 β / ( π M ) 0 . (6.11)For the dissipation functional, we obtain F ( Q, ˙ Q ) = 2 M mbβ (cid:90) d y sech ( βy ) sech( y + q ) tanh( y + q ) [ − ˙ θ + ˙ py − p ˙ q + ˙Λ y − y Λ ˙ q ] , (6.12)which gives w M = 0 , w θ = − M mbβ f ( β, q ) , w q = − M mbβ [ p f ( β, q ) + 2Λ f ( β, q ) ] , (6.13) w p = 2 M mbβ f ( β, q ) , w β = 0 , w Λ = 2 M mbβ f ( β, q ) . For H ( Q ), using the 6-parameter Ansatz we now obtain H ( Q ) M = p + β π Λ β + 2 β pbm I ( β, q ) + 4 β bm Λ I ( β, q ) (6.14) − gM β − β (cid:104) b + m − − b m (cid:105) I ( β, q ) . All the integrals are defined in Appendix A. For v µ ( Q ) = ∂ µ H ( Q ) we obtain v M = p + β π Λ β + 2 β pbm I ( β, q ) + 4 β bm Λ I ( β, q ) (6.15a) − gM β − (cid:104) b + m − − b m (cid:105) β I ( β, q ) ,v θ = 0 , (6.15b) v q = − β M pbm f ( β, q ) − β M bm Λ f ( β, q ) (6.15c)+ M (cid:104) b + m − − b m (cid:105) β f ( β, q ) ,v p = 2 M p + 2
M β bm I ( β, q ) , (6.15d) v β = 2 M β − M π Λ β + 2 M pbm I ( β, q ) + 4 M bm Λ I ( β, q ) (6.15e) − gM − M (cid:104) b + m − − b m (cid:105) I ( β, q ) − M β pbm f ( β, q ) − M β bm Λ f ( β, q ) + (cid:104) b + m − − b m (cid:105) M β f ( β, q ) ,v Λ = 2 π M Λ3 β + 4 βM bm I ( β, q ) . (6.15f)The symplectic equations of motion are ˙ Q µ = f µν ( Q ) u ν ( Q ) , (6.16)from which we find: ˙ M = − M mbβf ( β, q ) , (6.17a)˙ θ = − p + 23 β − gβM + 13 mbpβ I ( β, q ) + 2 mbβ Λ I ( β, q ) (6.17b)+ 2 mbpβ f ( β, q ) − mbpβ f ( β, q ) − mbβ Λ f ( β, q ) − (cid:104) b + m − − b m (cid:105) β [ 3 I ( β, q ) − β f ( β, q ) ]˙ q = 2 p + 2 β mb I ( β, q ) − mbβ f ( β, q ) , (6.17c)˙ p = 23 mbβ p f ( β, q ) − mbβ Λ f ( β, q ) − (cid:104) b + m − − m b (cid:105) β f ( β, q ) (6.17d)˙ β = − mb β f ( β, q ) − β Λ − β π mb I ( β, q ) + 12 β π mb f ( β, q ) , (6.17e)˙Λ = − + 4 β π + 4 π β pmb I ( β, q ) + 8 π β Λ mb I ( β, q ) (6.17f) − gβ Mπ − β π (cid:104) b + m − − b m (cid:105) f ( β, q ) − β π bm p f ( β, q ) − β π bm Λ f ( β, q ) . In Eq. (6.17f), we use the identity (A10). Here M ( t ) is adynamic variable. In order for the variational trial wavefunction to match the exact solution at t = 0, the initial conditions are: q = 0 , p = 0 , β = 1 , Λ = 0 , θ = − t (6.18) gM = (4 b − m − . As a check, the right-hand-sides of Eqs. (6.17) vanish[except for ˙ θ (0) = −
1] at these initial values, which guar-antees that the exact solution is stationary. For non-zerovalues of q and/or β , the values of p and Λ are some-times fixed by setting ˙ q = 0 and ˙ β = 0, and solvingEqs. (6.17c), and (6.17e) for p and Λ , which gives: p = 12 (cid:104) ˙ q − mbβ I ( β , q ) (6.19a)+ 2 mbβ f ( β , q ) (cid:105) , Λ = 14 β (cid:104) − ˙ β − mbβ f ( β , q ) (6.19b) − π mbβ I ( β , q ) + 12 π mbβ f ( β , q ) (cid:105) . When m = 0, the external potential is real and ˙ M = 0.The stability of the solutions to this equation for arbi-trary κ and for repulsive and attractive potential V aswell as positive and negative g was studied using a va-riety of methods, and the stability properties and smalloscillation frequencies for q, p, β, Λ were determined inRef. [38]. For that problem when we set κ = 1 and m = 0, our equations simplify to˙ q = 2 p , (6.20)˙ β = − β Λ , ˙Λ = − + 4 β π − gβ Mπ − β π (cid:104) b − (cid:105) f ( β, q )which agrees with the results in Ref. [38] once we use thefact that f [ G, q, γ ] in that paper is just β f ( β, q ) here.At m = 0 the small oscillation equations for β and q decouple. Using the expansions of the integrals found inAppendix A, we find that the small oscillation equationsare: δ ˙ q = 2 δp (6.21) δ ˙ p = −
815 ( b − / δq so that δ ¨ q + ω q δq = 0 , (6.22) ω q = 1615 ( b − / . This agrees with the result from the 4-parameter Ansatz.However, we get a different frequency for the β oscilla-tion, δ ˙ β = − δ Λ , (6.23) δ ˙Λ = (cid:104) b
15 + 4 π − (cid:105) δβ , so that δ ¨ β + ω β δβ = 0 , (6.24) ω β = 4 (cid:104) b
15 + 4 π − (cid:105) . Plots of ω q and ω β for m = 0 are shown in Fig. 2(a). VII. LINEAR RESPONSE RESULTS FOR THESIX CC APPROXIMATION
We linearize the set of equations given in (6.17) byexpanding the equations about the exact solutions, Q µ = Q µ + δQ µ keeping only the first order terms. Note that Q µ are given in Eqs. (6.18). Using the expansions ofAppendix ?? , we find δ ˙ M = − π mb M δq , (7.1a) δ ˙ θ = − gδM + 7 π mb δp (7.1b)+ 13 (cid:104) π − (cid:16)
12 + π (cid:17) gM (cid:105) δβ ,δ ˙ q = π (cid:0) π − (cid:1) mb δβ + 2 δp , (7.1c) δ ˙ p = 415 [ gM − δq − π mb δ Λ , (7.1d) δ ˙ β = (cid:104) π − π (cid:105) mb δq − δ Λ , (7.1e) δ ˙Λ = 2 bm π δp − gπ δM (7.1f)+ 215 (cid:104) − gM + 30 π + 4 (cid:105) δβ , where we have used the relation, b + m − b m −
14 = 2 − gM . (7.2)Equations (7.1) are written as δ ˙ Q µ = M µν ( Q ) δQ ν , (7.3)from which we find: δ ¨ Q µ + W µν ( Q ) δQ ν = 0 (7.4) W µν ( Q ) = − M µσ ( Q ) M σν ( Q ) . Here W µν ( Q ) is Hermitian. The square of the linearizedoscillation frequencies ω are given by the eigenvalues of W µν ( Q ). One can show that the matrix W µν ( Q ) canbe split into two blocks, one of them coupling ( δq, δ Λ , δθ ),the other coupling ( δp, δβ, δM ). Both of these blocks giveidentical eigenvalues, a zero eigenvalue and two non-zeroeigenvalues. For example, using Eqs. (7.1), we find δ ¨ q − [ A δq + B δ
Λ ] = 0 , (7.5a) δ ¨Λ − [ D δq + E δ
Λ ] = 0 , (7.5b)0 (a) m = 0 (b) m = 1 . m = 1 . FIG. 2. Plots of the linear response frequencies ω q (lower curve in red) and ω β (upper curve in blue) as a function of b for (a) m = 0, (b) m = 1, and (c) m = 1 .
25. The black dots represent data from the numerical simulation (see Sec. IX). The product mb controls the strength of the imaginary part of the potential. where A = 815 ( gM −
4) (7.6a)+ (9 π − π − b m ,B = (cid:0) − π (cid:1) π bm , (7.6b) D = bm (cid:110) gM π + 2(3 π − π (7.6c)+ ( gM − − π )15 π (cid:111) ,E = − π − m b + 815 ( gM − , (7.6d)from which we find ω = 12 (cid:20) − ( A + E ) ± (cid:113) ( A − E ) + 4 BD (cid:21) . (7.7)Although these two frequencies increase together when m = 0 as a function of b , once we get near m = 1 they start repelling each other and the dependence of the lowerfrequency has a maximum as a function of b instead ofmonotonically increasing. This is shown in Fig. 2. Notethat when m = 0 and b < /
4, the potential becomes re-pulsive, which leads to ω < VIII. COMPUTATIONAL STRATEGY
In our previous sections we were able to develop a sixparameter variational approach to the time evolution ofslightly perturbed solutions of the NLS equation in an ex-ternal complex potential. We were able to get an explicitanalytic expression as a function of m, b , of two oscilla-tion frequencies that affect the response of the solutionto small perturbations. So the first question we wouldlike to answer is how does this analytic result compareto the actual response found by numerically solving theNLS equation. The second question we want to answer1is the domain of applicability of the variational approachin terms of predicting the actual time evolution of thelow order moments of the solution. This has two parts:(i) for fixed m, b how long is the approximation valid and(ii) as we increase the size of the complex part of thepotential, by say varying m for fixed b , when does thisapproach start losing its validity. In our approximationfor all b, m that correspond to an attractive potential,there is no translational instability. So we would like tosee in our numerical simulations, that for the case m = 1(and κ = 1), the translational instability that arises dueto mixing of the solution we are considering with the firstexcited state in the potential occurs at times much laterthan the domain of applicability of the six CC method.For that case when 0 . < b < .
37 a late time transla-tional instability was found.To study numerically the evolution of Eqs. (1.1), wehave used a homemade code using a Crank-Nicolsonscheme [40]. In Ref. [33] we have shown that the Crank-Nicolson scheme is a reliable method for successfully solv-ing Eq. (1.1) in the presence of a complex potential. Forthe sake of comparison with the analytical calculations,the initial soliton shape ψ ( x,
0) in our simulations is givenby Eqs. (6.1) and (6.2) at t = 0. The complex solitonshape in the transverse spatial domain x was representedin a regular grid with mesh size ∆ x = 2 × − andfree boundary conditions were imposed. The mesh sizewas chosen to be much smaller than the initial solitonwidth parameter 1 /β (0) = 1, so that very small vari-ations of the soliton position could be accurately mea-sured by using a center of mass definition, i.e. q = (cid:104) x (cid:105) ,where the expectation value is defined in Eq. (6.6). Thesoliton width W ( t ) is the square root of the normalizedsecond moment G = π / (12 β ( t )). The soliton widthparameter 1 /β ( t ) in the simulations was calculated byusing the expression 1 /β ( t ) = (cid:112) G ( t ) /G (0). The otherCC’s measured in the simulations were the amplitude A ( t ) = max x ∈ R (cid:112) ρ ( x, t ) and the mass M ( t ) given byEq. (3.1). IX. COMPARISON OF COLLECTIVEVARIABLE THEORIES WITH SIMULATIONS
Our potential is symmetric in b ↔ m . When either b or m = 0 the potential is real and the small oscilla-tion equations for q, p β, Λ decouple giving rise to sep-arate oscillation frequencies in that regime. Once theimaginary part turns on, we expect that these two oscil-lation frequencies appear to a certain degree in all thecollective coordinates. Note that in the collective coor-dinate approach the mass is related to the height andthe CC parameter β and is not an independent param-eter, namely M ( t ) = 2 A ( t ) /β ( t ). First let us choose g = 1 , κ = 1 , m = 0 and b = 1 to see how well our CC ap-proximation works when compared with numerical simu-lations when the potential is real. For our simulations wechoose the parameters g = 1 , κ = 1 , q = 0 . , β = 1 . T (cid:32) T (cid:32) T (cid:32) (a) Position q ( t ) vs. t T (cid:32) T (cid:32) (b) Width 1 /β ( t ) − /β (0) vs. t FIG. 3. Comparison of the four CC (blue line), six CC(red line) and numerical simulation (black line). Parametersand initial conditions are m = 0 , q = 0 . , β = 1 . M = M , we display only q and 1 /β − /β . For this choice thetwo linear response periods are T q = 7 . T β = 4 . and all other parameters those of the exact solution. Thesmall oscillation theory for this case predicts separateoscillation frequencies for q and β , namely T q = 7 . ω q = 0 .
800 and T β = 4 . ω β = 2 . /β = 1 from 1 /β to show the oscillationin the numerical simulations. We see that for q ( t ), boththe amplitude as well as period of oscillation are well re-produced by the six CC theory. This is shown in Fig. 3.For the width parameter 1 /β , the oscillation periodis 4 .
00 which agrees well with the linear response result4 . . .
55. Moreover,the soliton amplitude A ( t ) has the period 3 .
85, which israther close to the above value of 4 . g = 1 , κ = 1 , m = 1, and three values for2 b . First we choose b = 0 . b = 0 . < b < .
56; and b = 1 .
45 is located in the upperstability regime 1 . < b < . q = 0 , p = 0 , β = 1 , Λ = 0 , θ = 0, and gM = (4 b − m − /
18, see Eqs. (6.18). In orderto test the stability of the exact solution, we choose ICsthat are slightly different from the above values. Thisexcites intrinsic oscillations of the soliton which are seenin the time evolution of the CCs, which is obtained bysolving the six CC equations, Eqs. (6.17), by a Math-ematica program. These oscillations are compared withthe oscillations which are observed in the simulations, i.e.in the numerical solution of the NLS equation. In par-ticular, the frequencies, periods, and amplitudes of theoscillations are compared.For the case b = 0 . T CC = T CC = 7 .
14, compared to T sim = 7 .
69. Thismeans that the error in the CC theories is only 7%.For the case b = 0 . T CC = 5 . T CC = 6 .
25 and T sim = 6 .
67, the error is 6%.For the case b = 1 .
45 the four CC result poorly fitsthe numerical result. The six CC result is very anhar-monic and the oscillation amplitudes do not agree wellwith the simulations as seen in Fig. 6. Nevertheless, theperiods T CC = 8 .
33 and T sim = 7 .
69 agree within anerror of 8%. Interestingly, the spectra exhibit a secondfrequency which is obtained also in the linear responsetheory. Fig. 2(b) shows the two frequencies for all valuesof b . However, the simulations show only one frequency.So far we have always taken q = 0 . p gives very similar results, because the q and p oscillationsare related, see the relations below Eq. (6.21). Let usnow consider b = 1 .
45 and finite values for Λ whichwill also affect the width 1 /β because their oscillationsare related. Choosing a very small, negative value Λ = − . = − . = +0 . M ( t ) exhibits a secondpeak at T = 2 .
38 which is stronger than the first peak at T = 8 .
33. This second peak belongs to the upper branchin Fig. 2(b) which was obtained by our linear responsetheory. However, this peak is not seen in the simulations. (a) Position q ( t ) vs. t (b) Amplitde A ( t ) vs. t (c) Width 1 /β ( t ) vs. t FIG. 4. Comparison of the four CC (blue line), six CC (redline), and numerical simulation (black line) for b = 0 .
1. Pa-rameters and initial conditions are m = 1 , q = . q ( t ), A ( t ), and 1 /β ( t ). For this choice the two linearresponse periods are: T q = 7 .
025 and T β = 4 . X. CONCLUSIONS
In this paper we investigated the domain of applica-bility of a four and six collective coordinate approxima-tion to study the response of the nodeless solution ofthe NLS equation in the presence of a complex potentialto small perturbations. This type of approximation hadbeen used in the past to study the response of exact solu-3 (a) q ( t ) vs. t (b) A ( t )(c) 1 /β ( t ) FIG. 5. Comparison of the four CC (blue line), six CC (redline), and numerical simulation (black line) for b = 0 .
5. Otherparameters and initial conditions being the same as in Fig. 4. tions of the NLS equation when in the presence of weakharmonic complex external potentials. In this paper weinstead considered a PT -symmetric potential where wecould vary the strength of the complex part of the po- tential from zero to its maximum allowed value. Using asmall oscillation approximation to the CC equations wewere able to obtain analytic expressions for the two fre-quencies of small oscillation found in our six CC approx-imation. These frequencies were quite close to those thatwere found in the numerical simulations of the discretizedPDEs when we perturbed the initial conditions of the ex-act solution. This was true for all allowed values of theparameter product bm which governed the strength ofthe imaginary part of the potential. We found that as weincreased bm , the four CC approximation quickly brokedown. The six CC approximation was quite a reasonableapproximation even at bm = 1 /
2, but at the maximumvalue we studied bm = 1 .
45, it tracked accurately theposition of the solitary wave for less than 1/4 of a periodand then began to differ from the numerical solution.
ACKNOWLEDGMENTS
F.C. would like to thank the Santa Fe Institute andthe Center for Nonlinear Studies at Los Alamos NationalLaboratory for their hospitality. F.G.M. and N.R.Q.acknowledge financial support from the Ministerio deEconom´ıa y Competitividad (Spain) through FIS2014-54497-P. F.G.M. also acknowledges financial supportfrom the Plan Propio of Universidad de Seville and isgrateful for the hospitality of the Mathematical Instituteof the University of Seville (IMUS) and of the Theoret-ical Division and Center for Nonlinear Studies at LosAlamos National Laboratory. N.R.Q. also acknowledgesfinancial support from the Junta de Andalucia (Spain)under Projects No. FQM207 and the Excellent GrantP11-FQM-7276. E.A. gratefully acknowledges supportfrom the Fondo Nacional de Desarrollo Cientifico y tec-nologico (FONDECYT) project No. 1141223 and fromthe Programa Iniciativa Cientfica Milenio (ICM) GrantNo. 130001. A.K. is grateful to Indian National ScienceAcademy (INSA) for awarding him INSA Senior Scien-tist position at Savitribai Phule Pune University, Pune,India. B.M. and J.F.D. would like to thank the Santa FeInstitute for their hospitality. B.M. acknowledges sup-port from the National Science Foundation through itsemployee IR/D program. The work of A.S. was sup-ported by the U.S. Department of Energy.
Appendix A: Definition of integrals
We note that dd z sech ( z ) = − ( z ) tanh( z ) = − ( z ) sinh( z ) , (A1a)dd z tanh( z ) = sech ( z ) . (A1b)4 (a) q ( t ) vs. t (b) A ( t )(c) 1 /β ( t ) FIG. 6. Comparison of the four CC (blue line), six CC (red line), and numerical simulation (black line) for b = 1 .
45. Otherparameters and initial conditions being the same as in Fig. 4.
Some useful integrals are the following: (cid:90) d z sech ( z ) = 2 , (A2a) (cid:90) d z sech ( z ) = π , (A2b) (cid:90) d z sech ( z ) = 43 , (A2c) (cid:90) d z z sech ( z ) = π , (A2d) (cid:90) d z sech ( z ) tanh ( z ) = 23 . (A2e)We define: I ( β, q ) = (cid:90) d x sech ( βy ) sech( x ) = (cid:90) d y sech ( βy ) sech( y + q ) , (A3a) I ( β, q ) = (cid:90) d x sech ( βy ) sech ( x ) = (cid:90) d y sech ( βy ) sech ( y + q ) , (A3b) I ( β, q ) = (cid:90) d x y sech ( βy ) sech( x ) = (cid:90) d y y sech ( βy ) sech( y + q ) . (A3c)5Also, we define: f ( β, q ) = (cid:90) d y sech ( βy ) sech( y + q ) tanh( y + q ) , (A4a) f ( β, q ) = (cid:90) d y y sech ( βy ) sech( y + q ) tanh( y + q ) , (A4b) f ( β, q ) = (cid:90) d y y sech ( βy ) sech( y + q ) tanh( y + q ) , (A4c) f ( β, q ) = (cid:90) d y sech ( βy ) sech( y + q ) tanh( y + q ) , (A4d) f ( β, q ) = (cid:90) d y y sech ( βy ) sech( y + q ) tanh( y + q ) , (A4e) f ( β, q ) = (cid:90) d y sech ( βy ) sech ( y + q ) tanh( y + q ) , (A4f) f ( β, q ) = (cid:90) d y y sech ( βy ) tanh( βy ) sech ( y + q ) , (A4g) f ( β, q ) = (cid:90) d y y sech ( βy ) sech ( y + q ) tanh( y + q ) , (A4h) f ( β, q ) = (cid:90) d y y sech ( βy ) tanh( βy ) sech( y + q ) , (A4i) f ( β, q ) = (cid:90) d y y sech ( βy ) tanh( βy ) sech( y + q ) . (A4j)Partial derivatives of I ( β, q ) are given by ∂I ( β, q ) ∂q = − (cid:90) d y sech ( βy ) sech( y + q ) tanh( y + q ) = − f ( β, q ) , (A5a) ∂I ( β, q ) ∂β = − (cid:90) d y y sech ( βy ) tanh( βy ) sech( y + q ) = − f ( β, q ) . (A5b)Partial derivatives of I ( β, q ) are given by ∂I ( β, q ) ∂q = − (cid:90) d y sech ( βy ) sech ( y + q ) tanh( y + q ) = − f ( β, q ) , (A6a) ∂I ( β, q ) ∂β = − (cid:90) d y y sech ( βy ) tanh( βy ) sech ( y + q ) = − f ( β, q ) . (A6b)Partial derivatives of I ( β, q ) are given by ∂I ( β, q ) ∂q = − (cid:90) d y y sech ( βy ) sech( y + q ) tanh( y + q ) = − f ( β, q ) , (A7a) ∂I ( β, q ) ∂β = − (cid:90) d y y sech ( βy ) tanh( βy ) sech( y + q ) = − f ( β, q ) . (A7b)A useful identity is obtained by integration of f ( β, q ) by parts. Using ∂∂y sech ( βy ) = − β sech ( βy ) tanh( βy ) , (A8)we find − βf ( β, q ) = (cid:90) y sech ( y + q ) d (cid:8) sech ( βy ) (cid:9) (A9)= − (cid:90) sech ( βy ) d (cid:8) y sech ( y + q ) (cid:9) = − (cid:90) d y sech ( βy ) sech ( y + q ) + 2 (cid:90) d y sech ( βy ) sech ( y + q ) tanh( y + q )= − I ( β, q ) + 2 f ( β, q ) . I ( β, q ) − βf ( β, q ) = 2 f ( β, q ) . (A10)We use this identity in the ˙Λ equation, (6.17f). Next, we now consider the expansion of the integrals and find to firstorder: I (1 + δβ, δq ) = π − π δβ , (A11a) I (1 + δβ, δq ) = 43 − δβ , (A11b) I (1 + δβ, δq ) = − π δq , (A11c) f (1 + δβ, δq ) = π δq , (A11d) f (1 + δβ, δq ) = π π
48 (16 − π ) δβ , (A11e) f (1 + δβ, δq ) = π
48 ( −
32 + 3 π ) δq , (A11f) f (1 + δβ, δq ) = 815 δq , (A11g) f (1 + δβ, δq ) = 13 + 245 ( −
15 + π ) δβ , (A11h) f (1 + δβ, δq ) = 13 − π δβ , (A11i) f (1 + δβ, δq ) = π
96 (16 − π ) δq , (A11j) f (1 + δβ, δq ) = π π
32 ( −
16 + π ) δβ . (A11k) Appendix B: Generalized traveling wave method