Study of the stability of solitons for the non-linear Schrödinger equation

Study of the stability of solitons for the non-linear Schrödinger equation Fernando Carreño Navas Study of the stability of solitons for the non-linear Schrödinger equation Fernando Carreño Navas Memoria presentada como parte de los requisitos para la obtención del título de doble Grado en Física y Matemáticas por la Universidad de Sevilla. Tutorizada por Prof. Renato Álvarez-Nodarse Prof. Niurka Rodríguez Quintero Contents 1 Introduction 1 1.1 Wave equation ................................. 1 1.2 Real media ................................... 3 1.3 Dispersive media, wave trains and nonlinear equations .......... 4 1.4 Solitons ..................................... 5 1.5 Derivation of the generalized nonlinear Schrödinger equation ...... 6 2 Soliton of the NLS equation. Stationary solution 7 2.1 Conservation equation of probability flux .................. 7 2.2 Conditions for the stationary solution of the NLS equation ........ 7 2.3 Resolution of the Soliton equation ...................... 9 2.4 Invariance under Galilean transformations ................. 9 3 Stability 11 3.1 Eigenvalues and Eigenfunctions ....................... 12 3.1.1 Boundary Conditions ......................... 14 3.1.2 Convergence of the eigenfunctions at the origin .......... 16 3.2 Case 𝜆= 0 ................................... 17 4 Numerical resolution of the NLS equation 19 4.1 Stationary Soliton ............................... 19 4.2 Numerical analysis of stability ........................ 22 4.2.1 Numerical Analysis of 𝜅= 1∕4 ................... 23 4.2.2 Conclusions .............................. 25 5 Stationary solution of the parametrically driven, damped NLS equation 27 5.1 Conservation equation of probability flux .................. 27 5.2 Conditions for the stationary solution of the parametrically driven, damped NLS equation .................................. 28 6 Stability of the parametrically driven, damped NLS equation 31 6.1 Eigenvalues and eigenfunctions ....................... 32 6.1.1 Boundary conditions ......................... 34 7 Numerical stability analysis 37 7.1 Numerical Analysis .............................. 37 7.1.1 The soliton 𝜙−............................ 38 7.1.2 The soliton 𝜙+............................ 39 7.1.3 Conclusions .............................. 39 ii study of the stability of solitons for the non-linear schrödinger eqation 8 Conclusions 43 Abstract English Abstract The Nonlinear Schrödinger equation is a fundamental equation in various scientific disciplines, including optics, hydrodynamics, and information transmission using optical fibers. In this work, we study the existence and stability of solitons using analytical and numerical methods, presenting two distinct approaches and comparing their results. Initially, we explore soliton-type solutions in a modified Schrödinger equation, introducing a parameter, 𝜅, to modulate the nonlinearity of the medium and observe differences in the behaviour of solitons. Additionally, we introduce dissipation and a parametric force, which allow us to examine how solitons behave in a more realistic medium. The analysis of the problem reveals stability regions that depend on the new parameters corresponding to dissipation, a parametric force and nonlinearity. Resumen en Castellano La ecuación de Schrödinger no lineal es una ecuación esencial en varias disciplinas científicas, como son la óptica, hidrodinámica, así como la transmisión de información usando fibra óptica. En este trabajo, estudiamos la existencia y estabilidad de solitones utilizando métodos analíticos y numéricos, presentando dos enfoques distintos y comparando sus resultados. Inicialmente, estudiamos solitones en una ecuación de Schrödinger modificada, introduciendo un parámetro, 𝜅, para modular la no linealidad del medio y observar la diferencia en el comportamiento de los solitones. Además, introducimos disipación y una fuerza paramétrica, que nos permiten examinar cómo se comportan los solitones en un medio más realista. El análisis del problema revela regiones de estabilidad las cuales dependen de los nuevos parámetros correspondientes a la fuerza paramétrica y a la disipación, incluyendo, además la no linealidad. 1 Introduction Pure mathematics and physics are becoming ever more closely connected, though their methods remain different. One may describe the situation by saying that the mathematician plays a game in which he himself invents the rules while the physicist plays a game in which the rules are provided by Nature, but as time goes on it becomes increasingly evident that the rules which the mathematician finds interesting are the same as those which Nature has chosen. Paul Dirac1 The laws that govern our reality are usually defined by equations. For example, the Maxwell’s equations for Electromagnetism, Boltzmann’s equation for Statistical Mechanics, and Schrödinger and Dirac equations for the Quantum Physics, are very well known examples of that. From them, we could gain some insight into the world around us, but like many models used in Physics or in Mathematics, they are built under strong simplifications. Linear type equations give us an important tool to understand, in first order, how a complex system may behave, but it will not allow us to understand the insights of the system. In order to understand deeper how a system may behave, we need nonlinear equations, which are not as simple as linear equations, but they allow us to construct new models for reality. There are problems such as the study of water waves dynamics which cannot be explained using linear equations, or how electromagnetic waves propagate through a waveguide. In this work, we will study one of these nonlinear equations, the nonlinear Schrödinger equation, which is able to model how waves travel through a nonlinear medium. 1.1 Wave equation There are many examples of waves in our common experience: a perturbation of the surface of a lake produced by a stone, the music of your favorite band, or the light produced by your cellphone, are examples of waves. We are used to interact with many types of waves, and due to this, we do not think much about their origin. In fact, not so many people realize that the way we hear music or see a picture have so many similarities with an ocean wave. 1Lecture delivered on presentation of the JAMES SCOTT prize, February 6, 1939. Published in: Proceedings of the Royal Society (Edinburgh) Vol. 59, 1938-39, Part II pp. 122-129 8study of the stability of solitons for the non-linear schrödinger eqation where we have separated the complex function in one real positive function 𝜙that serves as the modulus and other real function Θwhich implements the angular part. We can rewrite equation (2.1) as follows: 𝑖(𝜙𝑡+ 2𝜙𝑥Θ𝑥+𝜙Θ𝑥𝑥)+(𝜙𝑥𝑥 −𝜙Θ𝑡−𝜙Θ2 𝑥+𝛾𝜙2𝜅+1)=0.(2.4) Due to 𝜙(𝑥, 𝑡)and Θ(𝑥, 𝑡)being real functions, we can separate the real and imaginary parts of the equation (2.4). Taking that into account, we get the following system of differential equations: 𝜙𝑡+ 2𝜙𝑥Θ𝑥+𝜙Θ𝑥𝑥 = 0, 𝜙𝑥𝑥 −𝜙Θ𝑡−𝜙Θ2 𝑥+𝛾𝜙2𝜅+1 = 0.(2.5) There are some conditions that must be taking into account to find the stationary solution •Constant probability density. 𝜕𝑁 𝜕𝑡 =𝜕|𝑢|2 𝜕𝑡 =𝜕𝜙2 𝜕𝑡 = 2𝜙𝜙𝑡= 0.(2.6) From this, we obtain that 𝜙does not depend on time. •No probability density flux. Because we are looking for a stationary solution, then there must be no flux =𝑖(𝑢⋆ 𝑥𝑢−𝑢𝑥𝑢⋆) = 2Θ𝑥𝜙2.(2.7) Thus Θ𝑥= 0, (𝜙(𝑥)≢0). •Finite energy. This corresponds to bounded solutions (and its derivative) which goes to zero as 𝑥→±∞: lim 𝑥→±∞ |𝑢|= lim 𝑥→±∞ 𝜙= 0,lim 𝑥→±∞ |𝑢𝑥|= lim 𝑥→±∞ 𝜙𝑥= 0.(2.8) Using the condition (2.6), (2.7) into (2.5), we obtain that 𝜙𝑥𝑥(𝑥) − 𝜙(𝑥)Θ𝑡(𝑡) + 𝛾𝜙2𝜅+1(𝑥)=0→ 𝜙𝑥𝑥(𝑥) + 𝛾𝜙2𝜅+1(𝑥) 𝜙(𝑥)= Θ𝑡(𝑡) = 𝜇. (2.9) Since we have a time-dependent function on one side and a position-dependent function on the other side, they must be equal to a real constant, denoted by 𝜇. From these equations, we can conclude that Θ(𝑡) = 𝜇𝑡 (we have chosen the constant to be zero, as it is a global phase in the wave function). Thus, we get the following equation for 𝜙: 𝜙𝑥𝑥(𝑥) + 𝛾𝜙2𝜅+1(𝑥) − 𝜇𝜙(𝑥)=0.(2.10) In conclusion, in order to obtain the stationary solution we need to solve the last ordinary differential equation (2.10) and its solutions must satisfy the conditions (2.6), (2.8). stationary soliton 9 2.3 Resolution of the Soliton equation We will proceed with the resolution of the differential equation by multiplying equation (2.10) with 𝜙𝑥, integrating the resulting equation over the real axis and using the boundary conditions (2.8). This yields 𝜙2 𝑥= − 𝛾 𝑘+ 1𝜙2𝜅+2 +𝜇𝜙2,(2.11) or, equivalently, to solve the integral ±∫𝑑𝜙 √−𝛾 𝜅+1 𝜙2𝜅+2 +𝜇𝜙2 =𝑥−𝐶, (2.12) where 𝐶is a constant, that can be taken as zero if we choose the appropriate reference system. The integral (2.12) results in −1 √𝜇𝜅 arcsech⎛⎜⎜⎝[(𝛾 (𝜅+ 1)𝜇)1∕2𝜅 𝜙]𝑘⎞⎟⎟⎠ = ±𝑥. From this, we find the following explicit expression for 𝜙, 𝜙=((𝜅+ 1)𝜇 𝛾)1∕2𝜅 sech1∕𝜅(𝜅√𝜇𝑥).(2.13) In order to simplify the expression (2.13), it is convenient to define the following constants: 𝜇=𝜂2 𝜅2, 𝐴 =((𝜅+ 1)𝜇 𝛾)1∕2𝜅 . Then, the solution of equation (2.1) takes the form 𝑢(𝑥, 𝑡) = 𝜙(𝑥)𝑒𝑖𝜇𝑡 =𝐴sech1∕𝜅[𝜂𝑥]𝑒𝑖𝜇𝑡.(2.14) 2.4 Invariance under Galilean transformations Finally, we will show that the stationary solution (2.14) is invariant under a Galileo transformation. This will allow us to speak of a traveling soliton, and thus send information using this way of communication. The equation (2.1) is invariant under Galileo transformations if the equation maintains the same structure when we realize the following change of variables: 𝑥=𝑥−𝑣𝑡, 𝑡 =𝑡, 𝑢 =𝑢𝑒−𝑖𝑣 2(𝑥−1 2𝑣𝑡).(2.15) 10 study of the stability of solitons for the non-linear schrödinger eqation Using (2.15), the derivatives take the following form (considering that 𝜕𝑥∕𝜕𝑥 = 1 and 𝜕𝑥∕𝜕𝑡 = −𝑣): 𝜕𝑢 𝜕𝑥 =𝑒 𝑖𝑣 2(𝑥−1 2𝑣𝑡)[𝑢𝑥+𝑖𝑣 2𝑢], 𝜕2𝑢 𝜕𝑥2=𝑒 𝑖𝑣 2(𝑥−1 2𝑣𝑡)[𝑢𝑥𝑥 +𝑖𝑣𝑢𝑥−(𝑣 2)2 𝑢], 𝜕𝑢 𝜕𝑡 =𝑒 𝑖𝑣 2(𝑥−1 2𝑣𝑡)[𝑢𝑡−𝑣𝑢𝑥−𝑖(𝑣 2)2 𝑢]. (2.16) Substituting (2.16) into equation (2.1), we get 𝑖𝑢𝑡+𝑢𝑥𝑥 +𝛾|𝑢|2𝜅𝑢= 0.(2.17) Therefore, equation (2.1) is invariant under a Galilean transformation. 3 Stability According to the previous definition of stability established at the end of the section 1.4, we will study how the soliton behaves when we introduce a perturbation. We consider a perturbation of the following form: Ψ(𝑥, 𝑡)=[𝜙(𝑥) + 𝜖𝐯(𝑥, 𝑡)]𝑒𝑖𝜇𝑡,(3.1) where 𝜖is a small parameter1,𝜙(𝑥)satisfies 𝜙𝑥𝑥 −𝜇𝜙 +𝛾𝜙2𝜅+1 = 0,(3.2) and 𝐯(𝑥, 𝑡)represents a bounded perturbation of the steady state, 𝜙. Taking into account that 𝜙(𝑥)is real, 𝐯(𝑥, 𝑡)satisfies 𝑖𝐯𝑡+𝐯𝑥𝑥 + [(𝜅+ 1)𝛾𝜙2𝜅−𝜇]𝐯+𝛾𝜅𝜙2𝜅𝐯⋆= 0.(3.3) Multiplying equation (3.3) by 𝐯⋆and adding its complex conjugate to the resulting expression, we obtain 𝜕 𝜕𝑡|𝐯|2+𝜕 𝜕𝑥 𝐯𝑥𝐯⋆−𝐯⋆ 𝑥𝐯 𝑖=𝑖𝛾𝜅𝜙2𝜅(𝑥)[(𝐯⋆)2− (𝐯)2],(3.4) which after integration become 𝜕 𝜕𝑡 ∫ℝ|𝐯|2𝑑𝑥 =𝑖𝛾𝜅 ∫ℝ 𝜙2𝜅(𝑥)[(𝐯⋆)2− (𝐯)2]𝑑𝑥, (3.5) whenever (𝐯𝑥𝐯⋆−𝐯⋆ 𝑥𝐯)|+∞ −∞ = 0. This means that, a priori the norm of 𝐯is not conserved. Defining 𝐯=𝑓(𝑥, 𝑡) + 𝑖𝑔(𝑥, 𝑡), where both 𝑓, 𝑔 are real functions, (3.3) transforms into 𝑓𝑥𝑥 +((2𝜅+ 1)𝛾𝜙2𝜅−𝜇)𝑓=𝑔𝑡, 𝑔𝑥𝑥 +(𝛾𝜙2𝜅−𝜇)𝑔= −𝑓𝑡,(3.6) and we can separate these functions as a product of a temporal and a spatial function, i.e., 𝑓(𝑥, 𝑡) = 𝑇𝑓(𝑡)𝑣(𝑥), 𝑔(𝑥, 𝑡) = 𝑇𝑔(𝑡)𝑢(𝑥),(3.7) and, therefore, 𝑣𝑥𝑥 +((2𝜅+ 1)𝛾𝜙2𝜅−𝜇)𝑣 𝑢=𝑇𝑔 𝑡 𝑇𝑓=𝛼1, −𝑢𝑥𝑥 +(𝛾𝜙2𝜅−𝜇)𝑢 𝑣=𝑇𝑓 𝑡 𝑇𝑔=𝛼2, (3.8) 1We require 𝜖to satisfy that ||𝜖𝐯||∞≪||𝜙||∞, where ||𝑤||∞∶= sup{|𝑤(𝑥, 𝑡)|, 𝑥 ∈ℝ, 𝑡 ∈ [0, 𝑇 ]} 12 study of the stability of solitons for the non-linear schrödinger eqation where 𝛼1, 𝛼2are real constants. The equations for the temporal part are 𝑇𝑔 𝑡=𝛼1𝑇𝑓, 𝑇 𝑓 𝑡=𝛼2𝑇𝑔,(3.9) which can be decoupled into 𝑇𝑔 𝑡𝑡 =𝛼1𝛼2𝑇𝑔, 𝑇 𝑓 𝑡𝑡 =𝛼1𝛼2𝑇𝑓.(3.10) The characteristic equation of (3.10) is 𝜆2−𝛼1𝛼2= 0 →𝜆±= ±√𝛼1𝛼2.(3.11) Therefore, the temporal behaviour will be determinate by the value of 𝛼1𝛼2, which will give us the value of the elements in the spectrum. From this, we see that 𝜆2must be real and so 𝜆must be real or pure imaginary. The general solutions of equations (3.10) are 𝑇𝑔(𝑡) = 𝐸𝑔𝑒𝜆𝑡 +𝐹𝑔𝑒−𝜆𝑡, 𝑇 𝑓(𝑡) = 𝐸𝑓𝑒𝜆𝑡 +𝐹𝑓𝑒−𝜆𝑡,(3.12) where 𝐸𝑔,𝐹𝑔,𝐸𝑓, and 𝐹𝑓are complex constants so that 𝑇𝑔(𝑡)and 𝑇𝑓(𝑡)are real functions. According to this analysis, the stationary solution is stable whenever ℜ(𝜆)≤0 and unstable if ℜ(𝜆)>0. From the spatial dependency we have that 𝑣𝑥𝑥 +𝑣((2𝜅+ 1)𝛾𝜙2𝜅−𝜇)=𝛼1𝑢, 𝑢𝑥𝑥 +𝑢(𝛾𝜙2𝜅−𝜇)= −𝛼2𝑣. (3.13) Let’s examine first the case for 𝛼1=𝛼2=𝜆= 0 𝑣𝑥𝑥 +𝑣((2𝜅+ 1)𝛾𝜙2𝜅−𝜇)= 0, 𝑢𝑥𝑥 +𝑢(𝛾𝜙2𝜅−𝜇)= 0.(3.14) Equations (3.14) have both analytical solutions, which are 𝜙𝑥and 𝜙respectively. 3.1 Eigenvalues and Eigenfunctions We can transform the system (3.13) using the following change of variables 𝑢=𝜙𝑢2, 𝑣 =𝜙𝑣2.(3.15) This results in the following system: 2𝜙′𝑢′ 2+𝜙𝑢′′ 2= −𝛼2𝜙𝑣2, 2𝜙′𝑣′ 2+𝜙𝑣′′ 2=𝛼1𝜙𝑢2− 2𝜅𝛾𝜙2𝜅+1𝑣2,(3.16) stability 13 that can be rearrange as [𝑑 𝑑𝑥 (𝜙2𝑑 𝑑𝑥)+ 2𝜅𝛾𝜙2𝜅+2 (0 0 0 1)](𝑢2 𝑣2)=𝜙2(0 −𝛼2 𝛼10)(𝑢2 𝑣2).(3.17) The reader must notice that the system (3.17) is invariant whether we change the variable 𝑥→−𝑥, which means that the eigenfunctions of the system can be chosen to be even or odd. We solve our problem in the region [0,+∞), with an appropriate boundary conditions at 𝑥= 0, which keep the parity of the solutions. It is not easy to work with the system (3.17), so we will make a change in the variable 𝑥,𝑦= sech(𝜂𝑥). (Note that the change of variables 𝑦= sech(𝜂𝑥)is injective in [0,+∞)). We get −𝑑 𝑑𝑦 (𝑦2 𝑘+1(1 − 𝑦2)1 2 𝑑𝑢2 𝑑𝑦 )=𝛼2 𝜂2𝑦2 𝑘−1(1 − 𝑦2)−1 2𝑣2, −𝑑 𝑑𝑦 (𝑦2 𝑘+1(1 − 𝑦2)1 2 𝑑𝑣2 𝑑𝑦 )= −𝛼1 𝜂2𝑦2 𝑘−1(1 − 𝑦2)−1 2𝑢2+ 2𝜅+ 1 𝜅𝑦2 𝑘+1(1 − 𝑦2)−1 2𝑣2. (3.18) These equations could be simplified as 𝑦(−(2 𝜅+ 2)𝑦2+(2 𝜅+ 1))𝑑𝑢2 𝑑𝑦 + (1 − 𝑦2)𝑦2𝑑2𝑢2 𝑑𝑦2= −𝛼2 𝜂2𝑣2, 𝑦(−(2 𝜅+ 2)𝑦2+(2 𝜅+ 1))𝑑𝑣2 𝑑𝑦 + (1 − 𝑦2)𝑦2𝑑2𝑣2 𝑑𝑦2= −2𝜅+ 1 𝜅𝑦2𝑣2+𝛼1 𝜂2𝑢2. (3.19) We will use the Frobenius method (see, e.g., [3, Sec 8.5 p. 454]) to expand the solutions of (3.19) 𝑢2=ℜ(𝑦𝛽−1 𝜅 ∞ ∑ 𝑛=0 𝐶𝜆 𝑛𝑦𝑛), 𝑣2=ℜ(𝑦𝛽−1 𝜅 ∞ ∑ 𝑛=0 𝐷𝜆 𝑛𝑦𝑛).(3.20) We will assume that 𝐶𝜆 0, 𝐷𝜆 0, are not both zero (if this is not the case we can redefine the constants and use other value of 𝛽). We obtain the following recurrence relationships using (3.20) in (3.19): 𝐴(𝑛)(𝐶𝜆 𝑛 𝐷𝜆 𝑛)=𝐵(𝑛)(𝐶𝜆 𝑛−2 𝐷𝜆 𝑛−2),(3.21) where 𝐴(𝑛) = ⎛⎜⎜⎜⎝ (𝑛+𝛽)2−1 𝜅2 𝛼2 𝜂2 −𝛼1 𝜂2(𝑛+𝛽)2−1 𝜅2⎞⎟⎟⎟⎠ ,(3.22) 𝐵(𝑛) = ⎛⎜⎜⎝ 𝑏(𝑛) + 𝜅+ 1 𝜅0 0𝑏(𝑛) − 𝜅+ 1 𝜅⎞⎟⎟⎠ ,(3.23) 14 study of the stability of solitons for the non-linear schrödinger eqation being 𝑏(𝑛) = (𝑛+𝛽−2− 1 𝜅)(𝑛+𝛽−1+ 1 𝜅)−𝜅+ 1 𝜅.(3.24) For solving the recurrence (3.21), the determinant of the matrix 𝐴(0) should be zero, because 𝐶𝜆 0, 𝐷𝜆 0are not both zero. This yields −𝛼1𝛼2= (𝜂2𝛽2−𝜇)2= −𝜆2.(3.25) From this, we can distinguish four different solutions according to their behaviour at infinity. We introduce the following quantities: 𝛽+,1= +1 𝜂√+𝑖𝜆 +𝜇, 𝛽+,2= +1 𝜂√−𝑖𝜆 +𝜇, 𝛽−,1= −1 𝜂√+𝑖𝜆 +𝜇, 𝛽−,2= −1 𝜂√−𝑖𝜆 +𝜇. (3.26) We need to impose some boundary conditions in order to be able to characterize the eigenfunctions (and the generalized eigenfunctions) of our problem, as we will see in the next section 3.1.1. Summarizing, in order to obtain a priori solutions (3.20) to the system (3.17), we need to impose some boundary conditions in order to be able to obtain the solutions for the problem (3.18). This yields to the following expansions in the interval [0,∞): 𝑢(𝑥) = ℜ(𝐴sech𝛽(𝜂𝑥) ∞ ∑ 𝑛=0 𝐶𝜆 𝑛sech𝑛(𝜂𝑥)), 𝑣(𝑥) = ℜ(𝐴sech𝛽(𝜂𝑥) ∞ ∑ 𝑛=0 𝐷𝜆 𝑛sech𝑛(𝜂𝑥)). (3.27) 3.1.1 Boundary Conditions There are two boundary conditions that our functions must satisfy B1. 𝑢(𝑥), 𝑣(𝑥)must be bounded when 𝑥→∞. This condition will imply that we must impose that ℜ(𝛽)≥0. B2. Functions with well-defined parity. In order to be able to make the substitution 𝑦= sech(𝜂𝑥), the eigenfunctions should have a certain parity. This will imply that at 𝑥= 0, we must impose that either 𝑢(𝑥) = 𝑣(𝑥) = 0, to have an odd solution or 𝑢𝑥(𝑥) = 𝑣𝑥(𝑥)=0, to have an even solution. According to the definition of 𝛽(3.26), we can distinguish the following cases: 𝛽+,1, 𝛽+,2, 𝛽−,1, 𝛽−,2have a non-zero real part; or 𝛽±,1(respectively 𝛽±,2) are pure imaginary or zero. stability 15 Continuous Spectrum Without loss of generality, we will suppose first that 𝛽±,1are pure imaginary or zero, which is equivalent to consider 𝑖𝜆 ∈ℝand also ℑ(𝜆)≥𝜇, where ℑdenotes the imaginary part. In this case 𝛽+,2>0and 𝛽−,2<0, and so we need to discard 𝛽−,2<0, because the solution produced is divergent. This implies that in this region (ℑ(𝜆)≥𝜇), there are three independent pair of functions (𝑢𝑗, 𝑣𝑗),𝑗= 1,2,3, respectively. In the analogous case, if 𝛽±,2is pure imaginary, we obtain the region ℑ(𝜆)≤−𝜇. The conditions of parity B2 must be satisfied by a linear combination of (𝑢𝑗, 𝑣𝑗),𝑗= 1,2,3, to have a solution of (3.20). In order to have an odd solution, 𝑢𝑜𝑑𝑑(𝑥) = 𝑎1𝑢1(𝑥) + 𝑎2𝑢2(𝑥) + 𝑎3𝑢3(𝑥)and 𝑣𝑜𝑑𝑑(𝑥) = 𝑎1𝑣1(𝑥) + 𝑎2𝑣2(𝑥) + 𝑎3𝑣3(𝑥), where 𝑎1, 𝑎2, 𝑎3are real constants, it is necessary for 𝑢𝑜𝑑𝑑(0) and 𝑣𝑜𝑑𝑑(0) to be zero, thus 𝑎1𝑢1(0) + 𝑎2𝑢2(0) + 𝑎3𝑢3(0) = 0, 𝑎1𝑣1(0) + 𝑎2𝑣2(0) + 𝑎3𝑣3(0) = 0.(3.28) In the case of an even solution, 𝑢𝑒𝑣𝑒𝑛(𝑥) = 𝑎1𝑢1(𝑥) + 𝑎2𝑢2(𝑥) + 𝑎3𝑢3(𝑥)and 𝑣𝑒𝑣𝑒𝑛(𝑥) = 𝑎1𝑣1(𝑥) + 𝑎2𝑣2(𝑥) + 𝑎3𝑣3(𝑥), it is necessary for 𝑢𝑒𝑣𝑒𝑛 𝑥(0) and 𝑣𝑒𝑣𝑒𝑛 𝑥(0) to be zero, which implies that 𝑎1𝑢1 𝑥(0) + 𝑎2𝑢2 𝑥(0) + 𝑎3𝑢3 𝑥(0) = 0, 𝑎1𝑣1 𝑥(0) + 𝑎2𝑣2 𝑥(0) + 𝑎3𝑣3 𝑥(0) = 0.(3.29) For the continuous spectrum, we could always choose non-zero 𝑎1, 𝑎2, 𝑎3which satisfy (3.28) or (3.29). This implies that the spectrum contains 𝜆such that 𝑖𝜆 ∈ℝand |𝜆|≥𝜇. Discrete Spectrum If 𝛽+,1, 𝛽+,2, 𝛽−,1, 𝛽−,2have a non-zero real part then only two of them produce bounded functions, according to their definition (3.26). Therefore, we can produce only two independent pair of solutions, (𝑢𝑗, 𝑣𝑗),𝑗= 1,2, that must satisfy the parity conditions B2. In order to have an odd solution, 𝑢𝑜𝑑𝑑(𝑥) = 𝑎1𝑢1(𝑥) + 𝑎2𝑢2(𝑥)and 𝑣𝑜𝑑𝑑(𝑥) = 𝑎1𝑣1(𝑥) + 𝑎2𝑣2(𝑥), where 𝑎1, 𝑎2are real constants, it is necessary for 𝑢𝑜𝑑𝑑(0) and 𝑣𝑜𝑑𝑑(0) to be zero, so 𝑎1𝑢1(0) + 𝑎2𝑢2(0) = 0, 𝑎1𝑣1(0) + 𝑎2𝑣2(0) = 0.(3.30) In the case of an even solution, 𝑢𝑒𝑣𝑒𝑛(𝑥) = 𝑎1𝑢1(𝑥) + 𝑎2𝑢2(𝑥)and 𝑣𝑒𝑣𝑒𝑛(𝑥) = 𝑎1𝑣1(𝑥) + 𝑎2𝑣2(𝑥), it is necessary for 𝑢𝑒𝑣𝑒𝑛 𝑥(0) and 𝑣𝑒𝑣𝑒𝑛 𝑥(0) to be zero, which implies that 𝑎1𝑢1 𝑥(0) + 𝑎2𝑢2 𝑥(0) = 0, 𝑎1𝑣1 𝑥(0) + 𝑎2𝑣2 𝑥(0) = 0.(3.31) 16 study of the stability of solitons for the non-linear schrödinger eqation The previous systems (3.30), (3.31) do not produce any non-zero function unless the determinant is zero. The values of 𝜆where one of the conditions (3.30) or (3.31) is satisfied lead to the discrete spectrum. In order to calculate numerically the discrete spectrum, we will approximate the values of 𝑢𝑗(0), 𝑣𝑗(0), 𝑢𝑗 𝑥(0), 𝑣𝑗 𝑥(0),𝑗= 1,2using a finite number of terms (𝑁𝑚𝑎𝑥) of the series (3.27). In the numerical sections 4.2 and 7.1, we will refer to this method as “mixed”. Moreover, we define a new pair of functions Λ𝑜𝑑𝑑(𝜆) = 𝑢1(0)𝑣2(0) − 𝑣1(0)𝑢2(0),Λ𝑒𝑣𝑒𝑛(𝜆) = 𝑢1 𝑥(0)𝑣2 𝑥(0) − 𝑣1 𝑥(0)𝑢2 𝑥(0).(3.32) These pair of functions are the determinants of the systems (3.30) and (3.31), respectively, and so the zeros of (3.32) give us the discrete spectrum. For 𝜅= 1, the only eigenvalue in the discrete spectrum is 𝜆= 0, which corresponds to 𝜙(𝑥), 𝜙𝑥(𝑥)solutions. For 𝜅= 1∕4, we are able to find three more points in the discrete spectrum (the zero eigenvalue is always in the discrete spectrum), these points are associated to the eigenvalue 𝑖𝜆 = {0.845} which is the only odd function and 𝑖𝜆 = {0.575,0.980}, which are associated to even functions. The graphics of these functions are found in sections 4.6, 4.7,4.8, where it is discussed. 3.1.2 Convergence of the eigenfunctions at the origin It is important to verify whether our functions (3.27) are indeed bounded. We must look at the behaviour of the coefficients 𝐶𝜆 𝑛and 𝐷𝜆 𝑛when 𝑛is sufficiently large, and 𝛼1, 𝛼2are negligible. At the limit where 𝑛→∞the two equations of the recurrence relations (3.21) becomes 𝐶𝜆 𝑛=⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ (𝑛− 2) (𝑛−1+ 2 𝜅) 𝑛(𝑛+2 𝜅)𝐶𝜆 𝑛−2,for even 𝑛, 0,for odd 𝑛. Therefore, the system is decoupled, and we are able to study the convergence of the series. Thus, we need to prove that the following sum converges ∑ 𝑛=0 𝐶𝜆 2𝑛.(3.33) We will use the Raabe criterion for convergence of series (see, e.g., [3, Sec 5.2 p. 286]), which establishes that the series, (3.33), is convergent if 𝜀= lim inf 𝑛→∞[𝑛(𝐶𝜆 2𝑛−2 𝐶𝜆 2𝑛 − 1)]>1. stability 17 Since, in our case, 𝜀= 3∕2, the series (3.33) is convergent, and therefore is the series obtained by equation (3.21). 3.2 Case 𝜆= 0 The special case in our study appears when the matrix of the system (3.21), has a double zero. This only occurs when 𝛽= 1∕𝜅,𝜆= 0. In this case, the system has a double degeneracy at this value, and so we can find two eigenfunctions which linear combination give us the general solution. These functions can easily be obtained from (3.21) using an initial, 𝐶0 0= 0, 𝐷0 0= 1, and 𝐶0 0= 1, 𝐷0 0= 0, where 𝐶, 𝐷 are constants. •𝐶0 0= 1, 𝐷0 0= 0. In this case we obtained that 𝑢2= 1 and 𝑣2= 0, and then 𝑢=𝜙(𝑥), 𝑣 = 0, which is known to satisfy the original equation. This solution is even. •𝐶0 0= 0, 𝐷0 0= 1. In this case we obtained that 𝑢2= 0 and that 𝑣2is such that 𝐷0 𝑛={(1 − 3 𝑛)𝐷0 𝑛−2,for even 𝑛, 0,for odd 𝑛. From the above recurrence relation for the coefficients of the 𝑣2, and after some straightforward but cumbersome calculations we found that 𝑣2satisfies the following ordinary differential equation (ODE): (1 − 𝑦2)𝑑𝑣2 𝑑𝑦 = −𝑦𝑣2. The solution of this ODE is 𝑣2=√1 − 𝑦2, thus 𝑣2(𝑥) = tanh(𝜂𝑥). Then, 𝑢= 0, 𝑣 =𝜕𝑥𝜙(𝑥), which is an odd function. 24 study of the stability of solitons for the non-linear schrödinger eqation -15.0 -7.5 0.0 7.5 15.0 x 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 u61+(x) Solution u(61+) u(61+) -15.0 -7.5 0.0 7.5 15.0 x 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 u61-(x) Solution u(61-) u(61-) -15.0 -7.5 0.0 7.5 15.0 x -2.00 -1.50 -1.00 -0.50 0.00 0.50 1.00 1.50 2.00 [u61+(x)-u61-(x)] (10-9) Difference u(61+)-u(61-) Figure 4.3: Eigenfunctions ℜ(𝑢)corresponding to 𝜆1+ (left), 𝜆1− (center) and the difference between them (right). Data: 𝐿= 50, 𝑁 = 501, 𝜇 = 1, 𝛾 = 1, 𝜅 = 1. -15.0 -7.5 0.0 7.5 15.0 x -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 v62+(x) Solution v(62+) v(62+) -15.0 -7.5 0.0 7.5 15.0 x -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 v62-(x) Solution v(62-) v(62-) -15.0 -7.5 0.0 7.5 15.0 x -4.0 -2.0 0.0 2.0 4.0 6.0 8.0 [v62+(x)+v62-(x)] (10-9) Sum v(62+)+v(62-) Figure 4.4: Eigenfunctions ℜ(𝑣)corresponding to 𝜆2+ left, 𝜆2− center and the sum of them right. Data: 𝐿= 50, 𝑁 = 501, 𝜇 = 1, 𝛾 = 1, 𝜅 = 1. Mixed Numerical 𝑖𝜆30.575 0.574 𝑖𝜆40.845 0.845 𝑖𝜆50.980 0.980 Table 4.1: Table of values for the discrete spectrum, 𝜅= 1∕4, with two different methods, mixed with 𝑁𝑚𝑎𝑥 = 105(second column) and numerical results (third column). There is a small difference with respect to the mixed method as we can see in table 4.1. The eigenfunctions associated to the eigenvalues 𝜆3, 𝜆4, 𝜆5, calculated with the recurrence expression (3.21), are represented in figures 4.6,4.7 and 4.8, respectively. Notice that both 𝑢and 𝑣(obtained from equations (3.21)) are real functions according to (3.7). numerical analysis 25 -5 -2.5 0 2.5 5 Real part #10-7 -1.5 -1 -0.5 0 0.5 1 1.5 Imaginary part Spectrum for 5 = 1/4 Continuous Spectrum Discrete Spectrum Figure 4.5: Spectrum for 𝜅= 1∕4. Data: 𝐿= 50, 𝑁 = 501, 𝜇 = 1, 𝛾 = 1. -40.0 -30.0 -20.0 -10.0 0.0 10.0 20.0 30.0 40.0 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 Solution u(63) u(63) -40.0 -30.0 -20.0 -10.0 0.0 10.0 20.0 30.0 40.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 Solution v(63) v(63) Figure 4.6: Eigenfunction 𝑢corresponding to 𝜆3(left), eigenfunction 𝑣corresponding to 𝜆3(right). Data: 𝐿= 50, 𝑁 = 501, 𝜇 = 1, 𝛾 = 1, 𝜅 = 1∕4. 4.2.2 Conclusions From the numerical and theoretical results, we have extracted some interesting conclusions regarding the stability of solitons in the NLS equation. Theoretically, and with numerical evidence, we have obtained that the continuous spectrum is located on the imaginary axis, in the region {𝑖𝜆 ∈ℝ,|𝜆|≥𝜇}. Concerning the discrete spectrum, we have found that it must be located on the imaginary or real axis of the complex plane, and that 𝜆= 0 is an eigenvalue with eigenfunctions 𝜙and 𝜙𝑥. We know that the eigenvalues (discrete spectrum) satisfies the system (3.28) or (3.29), which is solved with a numerical approach. 26 study of the stability of solitons for the non-linear schrödinger eqation -40.0 -30.0 -20.0 -10.0 0.0 10.0 20.0 30.0 40.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 Solution u(64) u(64) -40.0 -30.0 -20.0 -10.0 0.0 10.0 20.0 30.0 40.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 Solution v(64) v(64) Figure 4.7: Eigenfunction 𝑢corresponding to 𝜆4(left), eigenfunction 𝑣corresponding to 𝜆4(right). Data: 𝐿= 50, 𝑁 = 501, 𝜇 = 1, 𝛾 = 1, 𝜅 = 1∕4. -40.0 -30.0 -20.0 -10.0 0.0 10.0 20.0 30.0 40.0 -5.0 -4.0 -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 Solution u(65) u(65) -40.0 -30.0 -20.0 -10.0 0.0 10.0 20.0 30.0 40.0 -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 5.0 Solution v(65) v(65) Figure 4.8: Eigenfunction 𝑢corresponding to 𝜆5(left), eigenfunction 𝑣corresponding to 𝜆5(right). Data: 𝐿= 50, 𝑁 = 501, 𝜇 = 1, 𝛾 = 1, 𝜅 = 1∕4. The mixed method allows us to numerically verify that for 𝜅≤1, there are no eigenvalues on the real axis, and thus the soliton should always be stable. We also established the existence of a discrete spectrum in {𝑖𝜆 ∈ℝ,|𝑖𝜆|< 𝜇}for values of 𝜅≤1. Numerically, we find a complete agreement with the theoretical results. 5 Stationary solution of the parametrically driven, damped NLS equation Previously, we have studied how a soliton is formed and how it behaves when we have only a dispersive-nonlinear medium. Despite this, in real media we must deal with many more factors in order to have a soliton. The most important factor to consider is dissipation, which produces a loss of the energy of the wave (in this case our soliton) reducing its amplitude when it propagates. With this hypothesis, the soliton can not maintain its shape. To solve this problem, we need to add energy to the system to compensate the energyloss due to dissipation. This, with some extra hypothesis, will allow us to create a soliton that will be studied in this chapter. Here we will start with the NLS equation studied in section 3.1, and we will add some terms, namely, 𝑖𝑢𝑡+𝑢𝑥𝑥 +𝛾|𝑢|2𝜅𝑢=ℎ𝑒𝑖𝜔𝑡𝑢⋆−𝑖𝛿𝑢, (5.1) where ℎis the magnitude of the oscillatory force, 𝛿is the dissipative parameter and ⋆ denotes the complex conjugate of a number. 5.1 Conservation equation of probability flux From the previous equation (5.1), we can obtain the following conservation equation of probability flux: 𝜕|𝑢|2 𝜕𝑡 −𝑖𝜕 𝜕𝑥 (𝑢𝑥𝑢⋆−𝑢⋆ 𝑥𝑢)=𝑖ℎ cos(𝜔𝑡)(𝑢2−(𝑢⋆)2) +ℎsin(𝜔𝑡)((𝑢⋆)2+𝑢2)−2𝛿|𝑢|2.(5.2) In comparison with (2.2), we have a non-zero source term in the conservation equation. The source term will produce an energy gain or loss, which will depend on which effect is dominant, the oscillatory force or the dissipation. As in section 2.1, we define 𝑁=|𝑢|2and =𝑖(𝑢⋆ 𝑥𝑢−𝑢𝑥𝑢⋆), as the probability density and the probability density flux, respectively. We also define the source term as 𝑆=𝑖ℎ cos(𝜔𝑡)(𝑢2− (𝑢⋆)2) + ℎsin(𝜔𝑡)((𝑢⋆)2+𝑢2)−2𝛿|𝑢|2. 28 study of the stability of solitons for the non-linear schrödinger eqation Because we are looking for a stationary solution, the probability density may remain constant in time and there will be no flux of the probability density in the system. 5.2 Conditions for the stationary solution of the parametrically driven, damped NLS equation In order to solve the equation (5.1), we will use the following ansatz 𝑢(𝑥, 𝑡) = 𝜙(𝑥, 𝑡)𝑒𝑖Θ(𝑥,𝑡),(5.3) where, as before, we have separated the complex function in one real positive function 𝜙that serves as the modulus, and other real function Θwhich implements the angular part. We can rewrite the equation (5.1) as follows: 𝑖(𝜙𝑡+ 2𝜙𝑥Θ𝑥+𝜙Θ𝑥𝑥 −ℎ𝜙 sin (𝜔𝑡 − 2Θ)+𝛿𝜙) +(𝜙𝑥𝑥 −𝜙Θ𝑡−𝜙Θ2 𝑥+𝛾𝜙2𝜅+1 −ℎ𝜙 cos (𝜔𝑡 − 2Θ))= 0.(5.4) The fact that 𝜙(𝑥, 𝑡)and Θ(𝑥, 𝑡)are real functions, allow us to separate the real and imaginary parts of the equation (5.4), so we find 𝜙𝑡+ 2𝜙𝑥Θ𝑥+𝜙Θ𝑥𝑥 =ℎ𝐹 sin (𝜔𝑡 − 2Θ)−𝛿𝜙, 𝜙𝑥𝑥 −𝜙Θ𝑡−𝜙Θ2 𝑥+𝛾𝜙2𝜅+1 =ℎ𝜙 cos (𝜔𝑡 − 2Θ).(5.5) In this case, we must take into account the following conditions, in order to find the stationary solution: •Constant probability density. 𝜕𝑁 𝜕𝑡 =𝜕|𝑢|2 𝜕𝑡 =𝜕𝜙2 𝜕𝑡 = 2𝜙𝜙𝑡= 0.(5.6) From this, it follows that 𝜙does not depend on time. •No probability density flux. Because we are looking for a stationary solution, then there must be no flux =𝑖(𝑢⋆ 𝑥𝑢−𝑢𝑥𝑢⋆) = 2Θ𝑥𝜙2.(5.7) Thus, Θ𝑥= 0, (𝜙(𝑥)≢0). •Finite energy. This corresponds to bounded solutions (and its derivative) which goes to zero as 𝑥→±∞: lim 𝑥→±∞ |𝑢|= lim 𝑥→±∞ 𝜙= 0,lim 𝑥→±∞ |𝑢𝑥|= lim 𝑥→±∞ 𝜙𝑥= 0.(5.8) Substituting (5.6), (5.7) into (5.5), we get ℎ𝜙 sin (𝜔𝑡 − 2Θ)−𝛿𝜙 = 0, 𝜙𝑥𝑥 −𝜙Θ𝑡+𝛾𝜙2𝜅+1 =ℎ𝜙 cos (𝜔𝑡 − 2Θ).(5.9) stationary damped soliton 29 From the first equation in (5.9), we obtain a relationship between the oscillatory force and the dissipative term ℎsin (𝜔𝑡 − 2Θ)−𝛿= 0 →sin (𝜔𝑡 − 2Θ)=𝛿 ℎ.(5.10) From (5.10), we get two different solutions for Θ, Θ+=1 2(𝜔𝑡 − arcsin (𝛿 ℎ))=1 2(𝜔𝑡 −𝜃+), Θ−=1 2(𝜔𝑡 −𝜋+ arcsin (𝛿 ℎ))=1 2(𝜔𝑡 −𝜃−). (5.11) where 𝜃+= arcsin (𝛿∕ℎ)and 𝜃−=𝜋− arcsin (𝛿∕ℎ). According to the equation (5.11), we also get a relation between the parameters 𝛿and ℎ. They must satisfy that 𝛿≤ℎ. Substituting (5.11) in the second equation of (5.9), we find 𝜙𝑥𝑥(𝑥) − 𝜙(𝑥)Θ𝑡(𝑡) + 𝛾𝜙2𝜅+1(𝑥) = ℎ𝜙 cos (𝜔𝑡 − 2Θ), 𝜙𝑥𝑥(𝑥) − 𝜙(𝑥)𝜔 2+𝛾𝜙2𝜅+1(𝑥)=±ℎ𝜙√1 − (𝛿 ℎ)2 , (5.12) which can be simplified into 𝜙𝑥𝑥(𝑥) + 𝛾𝜙2𝜅+1(𝑥) − 𝜇±𝜙= 0,(5.13) where 𝜇±(which must be always positive) is given by 𝜇±=𝜔 2±√ℎ2−𝛿2=𝜔 2+ℎcos(𝜃±).(5.14) Equation (5.13) could be solved in the same way as equation (2.10). Therefore 𝑢±=𝐴±sech1∕𝜅(𝜂𝑥)𝑒𝑖Θ±(𝑡)=𝜙(𝑥)𝑒𝑖Θ±(𝑡),(5.15) 𝐴±=(𝜇±(𝜅+ 1) 𝛾)1∕2𝜅 , 𝜇±=𝜂2 ± 𝜅2.(5.16) Notice that the parametrically driven, damped NLS equation (5.1) is not invariant under a Galilean transformation, due to the oscillation force term. 6 Stability of the parametrically driven, damped NLS equation In this chapter, we will study the stability of the soliton-type solution (5.15) of equation (5.1). According to the definition of stability, we will consider a perturbation of the soliton with the following form: Ψ(𝑥, 𝑡)=[𝜙(𝑥) + 𝜖𝐯(𝑥, 𝑡)]𝑒𝑖Θ±(𝑡),(6.1) where 𝜖is a small parameter1,𝜙(𝑥)satisfies 𝜙𝑥𝑥 −𝜇±𝜙+𝛾𝜙2𝜅+1 = 0,(6.2) and 𝐯(𝑥, 𝑡)represents a bounded perturbation of the exact solution. Taking into account that 𝜙(𝑥)is real, 𝐯(𝑥, 𝑡)satisfies 𝑖𝐯𝑡+𝐯𝑥𝑥 + [(𝜅+ 1)𝛾𝜙2𝜅−𝜔∕2 + 𝑖𝛿]𝐯+𝛾𝜅𝜙2𝜅𝐯⋆=ℎ𝑒𝑖𝜃±𝐯⋆, −𝑖𝐯⋆ 𝑡+𝐯⋆ 𝑥𝑥 + [(𝜅+ 1)𝛾𝜙2𝜅−𝜔∕2 − 𝑖𝛿]𝐯⋆+𝛾𝜅𝜙2𝜅𝐯=ℎ𝑒−𝑖𝜃±𝐯.(6.3) Defining 𝐯=𝑓(𝑥, 𝑡) + 𝑖𝑔(𝑥, 𝑡), where both 𝑓, 𝑔 are real functions. (6.3) transforms into 𝑓𝑥𝑥 +((2𝜅+ 1)𝛾𝜙2𝜅−𝜇±)𝑓= 2𝛿𝑔 +𝑔𝑡, 𝑔𝑥𝑥 +(𝛾𝜙2𝜅−𝜔∕2 + ℎcos(𝜃±))𝑔= −𝑓𝑡.(6.4) As before, we will use the separation of variables 𝑓(𝑥, 𝑡) = 𝑇𝑓(𝑡)𝑣(𝑥), 𝑔(𝑥, 𝑡) = 𝑇𝑔(𝑡)𝑢(𝑥).(6.5) Then equation (6.3) transforms into 𝑣𝑥𝑥 +((2𝜅+ 1)𝛾𝜙2𝜅−𝜇±)𝑣 𝑢= 2𝛿𝑇𝑔 𝑇𝑓+𝑇𝑔 𝑡 𝑇𝑓=𝛼1, −𝑢𝑥𝑥 +(𝛾𝜙2𝜅−𝜔∕2 + ℎcos(𝜃±))𝑢 𝑣=𝑇𝑓 𝑡 𝑇𝑔=𝛼2, (6.6) where 𝛼1, 𝛼2are real constants. The equations for the temporal part are 𝑇𝑔 𝑡+ 2𝛿𝑇 𝑔=𝛼1𝑇𝑓, 𝑇 𝑓 𝑡=𝛼2𝑇𝑔,(6.7) 1We require 𝜖to satisfy that ||𝜖𝐯||∞≪||𝜙||∞, where ||𝑤||∞∶= sup{|𝑤(𝑥, 𝑡)|, 𝑥 ∈ℝ, 𝑡 ∈ [0, 𝑇 ]} 32 study of the stability of solitons for the non-linear schrödinger eqation or, equivalently, 𝑇𝑔 𝑡𝑡 + 2𝛿𝑇 𝑔 𝑡=𝛼1𝛼2𝑇𝑔, 𝑇 𝑓 𝑡=𝛼2𝑇𝑔.(6.8) From equations (6.8), we deduce that both 𝑇𝑓and 𝑇𝑔have the same temporal behaviour. The characteristic equation of 𝑇𝑔is 𝜆2+ 2𝛿𝜆 −𝛼1𝛼2= 0 →𝜆±= −𝛿±√𝛿2+𝛼1𝛼2,(6.9) so the temporal behaviour will be determined by the value of 𝛼1𝛼2which will give us the value of the elements in the spectrum. Consequently, it follows that (𝜆+𝛿)2must be real. We need to consider two different cases 𝜆= −𝛿+𝜆𝑅,𝜆𝑅∈ℝ, which could produce instability in the system if 𝜆𝑅> 𝛿 and 𝜆= −𝛿+𝑖𝜆𝐼,𝜆𝐼∈ℝ, which will not produce any instability due to the fact that 𝛿≥0. For the spatial dependency, we have 𝑣𝑥𝑥 +𝑣((2𝜅+ 1)𝛾𝜙2𝜅−𝜇±)=𝛼1𝑢, 𝑢𝑥𝑥 +𝑢(𝛾𝜙2𝜅−𝜔∕2 + ℎcos(𝜃±))= −𝛼2𝑣. (6.10) 6.1 Eigenvalues and eigenfunctions Now we use the same change of variables as in section 3.1,𝑢=𝑢2𝜙,𝑣=𝑣2𝜙, where 𝜙=𝐴sech1∕𝜅(𝜂𝑥). This leads to the following Sturm-Liouville problem: −𝛼1𝜙𝑢2+𝑣2𝑥𝑥𝜙+ 2𝜙𝑥𝑣2𝑥+ 2𝜅𝛾𝜙2𝜅+1𝑣2= 0,(6.11) 𝛼2𝑣2𝜙+𝑢2𝑥𝑥𝜙+ 2𝜙𝑥𝑢2𝑥+ 2𝜙ℎ cos(𝜃±)𝑢2= 0,(6.12) or, in a more compact notation [𝑑 𝑑𝑥 (𝜙2𝑑 𝑑𝑥)+ 2𝜙2(ℎcos(𝜃±) 0 0𝜅𝛾𝜙2𝜅)](𝑢2 𝑣2)=𝜙2(0 −𝛼2 𝛼10)(𝑢2 𝑣2).(6.13) As in section 3.1, we can restrict ourself to [0,∞) adding a proper boundary condition determined by the parity of the solution, which will be discussed in section 6.1.1. In the restricted dominion [0,∞), we can make the change of variables, 𝑦= sech(𝜂𝑥). The system (6.13) becomes 𝑦((1 + 2 𝜅)−(2 + 2 𝜅)𝑦2)𝑣2𝑦+𝑦2(1 − 𝑦2)𝑣2𝑦𝑦 + 2𝑣2 (𝜅+ 1) 𝜅𝑦2=𝛼1 𝑢2 𝜂2, 𝑦((1 + 2 𝜅)−(2 + 2 𝜅)𝑦2)𝑢2𝑦+𝑦2(1 − 𝑦2)𝑢2𝑦𝑦 + 2𝑢2 𝜂2ℎcos(𝜃±)=−𝛼2 𝑣2 𝜂2. (6.14) We will use the Frobenius method (see, e.g., [3, Sec 8.5 p. 454]) to expand the solutions of (6.14) 𝑢2=ℜ(𝑦𝛽−1 𝜅 ∞ ∑ 𝑛=0 𝐶𝜆 𝑛𝑦𝑛), 𝑣2=ℜ(𝑦𝛽−1 𝜅 ∞ ∑ 𝑛=0 𝐷𝜆 𝑛𝑦𝑛),(6.15) stability of the damping nls eqation 33 where 𝛽is a constant which depends on the value of 𝜆. With this expansion, we can write the following system of recurrence, which relates the coefficients 𝐶𝜆 𝑛and 𝐷𝜆 𝑛: 𝐴(𝑛)(𝐶𝜆 𝑛 𝐷𝜆 𝑛)=𝐵(𝑛)(𝐶𝜆 𝑛−2 𝐷𝜆 𝑛−2),(6.16) where now the matrices 𝐴(𝑛), 𝐵(𝑛)are given by 𝐴(𝑛) = ⎛⎜⎜⎜⎝ 𝜒(𝑛) + ℎ 𝜂2cos(𝜃±)𝛼2∕𝜂2 −𝛼1∕𝜂2𝜒(𝑛) − ℎ 𝜂2cos(𝜃±)⎞⎟⎟⎟⎠ ,(6.17) 𝐵(𝑛) = ⎛⎜⎜⎝ 𝑏(𝑛) + 𝜅+ 1 𝜅0 0𝑏(𝑛) − 𝜅+ 1 𝜅⎞⎟⎟⎠ ,(6.18) respectively, and 𝜒(𝑛)and 𝑏(𝑛)are defined by 𝜒(𝑛)=(𝑛+𝛽)2−1 𝜅2+ℎ 𝜂2cos(𝜃±), 𝑏(𝑛) = (𝑛+𝛽−2− 1 𝜅)(𝑛+𝛽−1+ 1 𝜅)−𝜅+ 1 𝜅. (6.19) In order to have non-zero solutions the determinant of 𝐴(0) should be zero, because of the condition that 𝐶𝜆 𝑛−2 =𝐷𝜆 𝑛−2 = 0. This implies that 𝜂4𝜒2(0) = ℎ2−𝛿2−𝛼1𝛼2→ℎ2−(−𝜂2𝛽2+𝜔 2)2 = (𝜆+𝛿)2.(6.20) From equation (6.20), we can obtain the following four possible values for 𝛽: 𝛽+,1= +1 𝜂√𝜔 2+√ℎ2− (𝜆+𝛿)2, 𝛽+,2= +1 𝜂√𝜔 2−√ℎ2− (𝜆+𝛿)2, 𝛽−,1= −1 𝜂√𝜔 2+√ℎ2− (𝜆+𝛿)2, 𝛽−,2= −1 𝜂√𝜔 2−√ℎ2− (𝜆+𝛿)2. (6.21) We need to use the boundary conditions B1 and B2 defined in section 3.1.1, to find the eigenfunctions (and the generalized eigenfunctions) of our problem (6.10). Summarizing, in order to obtain a priori solutions to our system (6.14), we need to impose some boundary conditions to obtain the solutions of problem (6.10). This yields the following expansions in [0,∞): 𝑢=ℜ(𝐴sech𝛽(𝜂𝑥) ∞ ∑ 𝑛=0 𝐶𝜆 𝑛sech𝑛(𝜂𝑥)), 𝑣=ℜ(𝐴sech𝛽(𝜂𝑥) ∞ ∑ 𝑛=0 𝐷𝜆 𝑛sech𝑛(𝜂𝑥)). (6.22) 40 study of the stability of solitons for the non-linear schrödinger eqation -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 Real part -5 -4 -3 -2 -1 0 1 2 3 4 5 Imaginary part Spectrum for 5 = 1, ! = 2, h = 0.5 / = 0.2 of ?+ Continuous Spectrum Discrete Spectrum Spurious values Figure 7.2: Numerical Spectrum obtained for the soliton 𝜙+. Data: 𝜅= 1,𝜔= 2,ℎ= 0.5, 𝛿= 0.2,𝛾= 1,𝑁= 601, 𝐿 = 100. Mixed Numerical 𝜆10.000 0.000 𝜆𝑐-0.237 -0.237 Table 7.3: Table of results for the mixed and numerical method for soliton 𝜙+. Data: 𝜅= 1,𝜔= 2,ℎ= 1.2,𝛿= 0.9,𝛾= 1. Mixed Numerical 𝜆𝑐0.163 0.163 Table 7.4: Table of results for the mixed and numerical method for soliton 𝜙+. Data: 𝜅= 1,𝜔= 2,ℎ= 1.2,𝛿= 0.5,𝛾= 1. Mixed Numerical 𝜆1-0.233 -0.232 𝜆2-0.064 -0.061 𝜆30.000 0.000 𝜆40.106 0.106 𝜆𝑐-0.237 -0.237 Table 7.5: Table of results for the mixed and numerical method for soliton 𝜙+. Data: 𝜅= 1∕4,𝜔= 2,ℎ= 1.2,𝛿= 0.9,𝛾= 1. the existence of two different solitons, 𝜙−and 𝜙+, each exhibiting a completely different behaviour. numerical stability analysis 41 -2 -1.5 -1 -0.5 0 Real part -5 -4 -3 -2 -1 0 1 2 3 4 5 Imaginary part Spectrum for 5 = 1, ! = 2, h = 1.2 / = 0.9 of ?+ Continuous Spectrum Discrete Spectrum Spurious values Figure 7.3: Numerical Spectrum obtained for the soliton 𝜙+. Data: 𝜅= 1,𝜔= 2,ℎ= 1.2, 𝛿= 0.9,𝛾= 1,𝑁= 601, 𝐿 = 100. -1.5 -1 -0.5 0 0.5 Real part -5 -4 -3 -2 -1 0 1 2 3 4 5 Imaginary part Spectrum for 5 = 1, ! = 2, h = 1.2 / = 0.5 of ?+ Continuous Spectrum Spurious values Figure 7.4: Numerical Spectrum obtained for the soliton 𝜙+. Data: 𝜅= 1,𝜔= 2,ℎ= 1.2, 𝛿= 0.5,𝛾= 1,𝑁= 601, 𝐿 = 100. The soliton 𝜙−is found to be unstable for every selection of parameters, because of the existence of a single eigenvalue with positive real part (see figure 7.1). The discrete eigenvalue could be calculated using the mixed or the numerical method to obtain approximate results. The continuous spectrum of this soliton is found analytically, which does not introduce any instability to the system, as the real part is always negative. 42 study of the stability of solitons for the non-linear schrödinger eqation -2 -1.5 -1 -0.5 0 Real part -5 -4 -3 -2 -1 0 1 2 3 4 5 Imaginary part Spectrum for 5 = 1/4, ! = 2, h = 1.2 / = 0.9 of ?+ Continuous Spectrum Discrete Spectrum Figure 7.5: Numerical Spectrum obtained for the soliton 𝜙+. Data: 𝜅= 1∕4,𝜔= 2, ℎ= 1.2,𝛿= 0.9,𝛾= 1,𝑁= 601, 𝐿 = 100. The soliton 𝜙+presents stability regions with a strong dependency in the parameters. The continuous spectrum could be analytically studied (see equation (6.23)), which can produce instability in the system for certain parameters set (see figure 7.4). Moreover, we find that there exists a discrete spectrum which can produce instability, similar to 𝜙− (see figure 7.5) or due to the appearance of spurious eigenvalues (see figure 7.3). Despite the regions of instability, there exists regions where the soliton 𝜙+is complete stable (see figure 7.2). Remark: The analytical results obtained do not support the existence of spurious eigenvalues. This discrepancy arises due to an incorrect ansatz, (6.5), which does not account for this case. Supported by a scholarship from IMUS (Institute of Mathematics of the University of Sevilla), we were able to develop a more general approach based on the ansatz proposed in [7]. The findings of this study have been published in [8], where numerical and analytical results for the two solitons are presented. For the soliton 𝜙+, we derived the complete ℎ-𝛿stability diagram for various values of 𝜅, ranging from 𝜅= 0.5 to 𝜅= 3.5. This provides a complete characterization of the stability of the soliton for each pair of ℎand 𝛿, extending the results of [7]. For the soliton 𝜙−, analytical results indicate that it is inherently unstable. 8 Conclusions The NLS equation is a fundamental equation that arises in many different fields such as optics, transmission lines, or water waves. Therefore, the results regarding the stability of its soliton-type solutions are valuable as they could significantly improve how information travels through a real medium. In this work, we analyze the existence and stability of solitons using both analytical and numerical approaches, presenting two different methods and comparing the obtained results. First, we analyze how to obtain soliton-type solutions for the modified Schrödinger equation, which includes a new parameter 𝜅, that allows us to introduce the nonlinearity of the medium and observe the emergence of different solitons. Next, we analytically solve the eigenvalue problem and found a recursive expression that allows us to construct the eigenfunctions. We also present a complete numerical treatment to find the spectrum and compare it with the analytical results. Using both numerical methods, we have obtained new evidence regarding the existence of a discrete spectrum of eigenvalues, whose values present a strong dependency on 𝜅. The used code has been developed in MATLAB, and for the sake of the space we have not include them here, nevertheless, they are available by request. We also consider the NLS equation with dissipation and an oscillatory force. As a novel result of this work, two stationary solutions are obtained and their stability is investigated. The obtained results are quite surprising due to the existence of different stability regions based on the parameters and the existence of two types of solitons 𝜙−and 𝜙+. These results will be the core of the publication [8] that is under preparation and that will be submitted to a Mathematical-Physics Journal. The problem of obtaining analytically the complete discrete spectrum remains as an open question, but we hope that this work helps to make a step further in the resolution of it. Bibliography [1] Alexeeva, N. V., Barashenkov, I. V., & Pelinovsky, D. E., Dynamics of the parametrically driven NLS solitons beyond the onset of the oscillatory instability. Nonlinearity, 12(1), (1999), 103-140. https://doi.org/10.1088/0951-7715/12/1/007 [2] Allaire, G., Numerical Analysis and Optimization: An Introduction to Mathematical Modelling and Numerical Simulation. Oxford University Press, Oxford, 2007. [3] Arfken, G. B., Mathematical Methods for Physicists. Academic Press, New York, 1985. [4] Balakrishnan, R., Soliton propagation in nonuniform media. Physical Review A, 32(2), (1985), 1144-1149. https://doi.org/10.1103/physreva.32.1144 [5] Balanis, C. A., Antenna theory: Analysis and Design. John Wiley & Sons, New Jersey, 2005. [6] Ballentine, L. E., Quantum Mechanics: A Modern Development (2nd Edition). World Scientific Publishing Company, Singapore, 2014. [7] Barashenkov, I. V., Bogdan, M. M., & Korobov, V. I., Stability Diagram of the Phase- Locked Solitons in the Parametrically Driven, Damped Nonlinear Schrödinger Equation. Europhysics Letters, 15(2), (1991), 113-118. https://doi.org/10.1209/ 0295-5075/15/2/001 [8] Carreño-Navas, F., Álvarez-Nodarse, R., & Quintero N. R., Oscillatory instability and stability of stationary solutions in the parametrically driven, damped nonlinear Schrödinger equation. https://doi.org/10.48550/arXiv.2411.17655. [9] R. Carretero-González, D.J. Frantzeskakis, and P.G. Kevrekidis. Nonlinear Waves & Hamiltonian Systems: From One To Many Degrees of Freedom, From Discrete To Continuum. Oxford University Press, 2024. [10] Griffiths, D. J. Introduction to Electrodynamics. Pearson Higher Education, Glenview, IL, 2014. [11] Hansen, P. J., & Nicholson, D. R. Simple soliton solutions. American Journal Of Physics, 47(9), (1979), 769-771. https://doi.org/10.1119/1.11925 [12] Knoll, G. F., Radiation Detection and Measurement. John Wiley & Sons, New Jersey, 2010. [13] Quintero, N. R., Mertens, F. G., & Bishop, A. R. (2015). Soliton stability criterion for generalized nonlinear Schrödinger equations. Physical Review E, 91(1), 012905. https://doi.org/10.1103/physreve.91.012905 46 study of the stability of solitons for the non-linear schrödinger eqation [14] Remoissenet, M. Waves called solitons: Concepts and Experiments. Springer Science & Business Media, New York, 2013. [15] Scott, A. Nonlinear science: Emergence and Dynamics of Coherent Structures. Oxford University Press, Oxford, 2003. [16] Vitanov, N. K., Chabchoub, A., & Hoffmann, N. Deep-Water Waves: on the Nonlinear Schrödinger Equation and its Solutions. Journal Of Theoretical And Applied Mechanics, 43(2), (2013), 43-54. [17] Wang P., & Li JP. Quantum Interpretation of Frequency Operator. In ICNC07: Third Internation Conference on Natural Computation. China: IEEE Computer Society Press, (2007), 613-618. [18] Yang, J., & Nixon, S., Stability of soliton families in nonlinear Schrödinger equations with non-parity-time-symmetric complex potentials. Physics Letters A, 380(45), 3803-3809. https://doi.org/10.1016/j.physleta.2016.09.023