Pseudo Spectral Transform for Schr\"odinger--Poisson Equations

TL;DR

This paper introduces a novel pseudo spectral transform method to find exact stationary solutions to the Schr"odinger--Poisson equations, addressing technical challenges through a generalized spectral approach and integral convergence analysis.

Contribution

It develops a new pseudo spectral transform technique for solving Schr"odinger--Poisson equations, extending the Korteweg--de Vries kernel to overcome previous methodological limitations.

Findings

Derived exact stationary solutions using the pseudo spectral transform.

Established convergence of the integral representation.

Compared the new method with homotopy analysis approach.

Abstract

Here we present exact, stationary, parametric solutions to the Schr\"odinger--Poisson system. We confront two images: on one hand, we draw on the homotopy analysis method which leads us to a nonlinear integral scheme. Indeed, this approach might be simplified by looking for sufficiently smooth solutions vanishing asymptotically. However, since our system possesses stiffness an additional analysis has to be considered. On the other hand, we seek for exact solutions over the inverse scattering method by introducing a pseudo spectral transform. In fact, this pseudo spectral method generalises Korteweg--de Vries family's kernel and let us to circumvent some technical difficulties originally arisen in our first approach although, again, we come to an integral representation, which we test for convergence.

Equations

$$\left[-△+\varphi \right]\psi$$

$$△\varphi$$

$$R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R=8\pi ⟨\psi |{\mathbf{T}}_{\mu \nu }|\psi ⟩,$$

$$\left[-{\partial }_x^2+b\varphi \right]u$$

$${\varphi }_{xx}$$

$$-u_{xxxx}+b(\varphi u_{xx}+{\varphi }_{xx}u)+2b{\varphi }_x{u}_x=\omega u_{xx},$$

$${\varphi }_x={\int }_{\mathcal{X}}\mathrm{d}x|u{|}^2=\frac{u_{xxxx}-u^{-1}(u_{xx}{)}^2-bu|u{|}^2}{2b{u}_x}.$$

$$u_{xxxx}=\frac{b|u{|}^2{u}^3+[(u_{xx}{)}^2+2u_x{u}_{xxx}]u-2(u_{xx}{)}^2{u}_x}{u^2},$$

$$\underset{x\to \pm \infty }{\lim }\left[|x{|}^{ϵ+1}u(x)\right]=0,0<ϵ,$$

$$u_{xxxx}\approx b{u}^3,$$

$$\begin{array}{rl}u(x)=C_1\frac{b{x}^3}{6}& {\int }^x\mathrm{d}y{u}^3(y)-C_2\frac{b{x}^2}{2}{\int }^x\mathrm{d}yy{u}^3(y)\\ +C_3\frac{bx}{2}& {\int }^x\mathrm{d}y{y}^2{u}^3(y)-C_4\frac{b}{6}{\int }^x\mathrm{d}y{y}^3{u}^3(y)+F(x),\end{array}$$

$${\Xi }_{\mu }(u^3)=\sum _{\nu =0}^{\mu }{u}_{\mu -\nu }\sum _{\beta =0}^{\nu }{u}_{\beta }{u}_{\nu -\beta },$$

$$u(x)=\sum _{\mu =0}^{\infty }{u}_{\mu }(x)=\sum _{\mu =0}^{\infty }\left\{{\int }^x\mathrm{d}yK(x,y){\Xi }_{\mu -1}(u^3)\right\},$$

$$\left|1-{\int }_{{\chi }_1}^{{\chi }_2}\mathrm{d}yK(y,y)\right|\le 1,$$

$$u(x)=\mathrm{sech}(ax)-b\sum _{\mu =1}^{\infty }\left\{\frac{C_2{x}^2}{2}{\int }^x\mathrm{d}yy{\Xi }_{\mu -1}(u^3)+\frac{C_4}{6}{\int }^x\mathrm{d}y{y}^3{\Xi }_{\mu -1}(u^3)\right\};$$

$$\mathbf{H}\Psi =k^2\Psi ,$$

$$\underset{x\to \infty }{\lim }[\exp (p_{\alpha }x){\Psi }_{\alpha }(x)]=1,\alpha =1,2,\dots ,A.$$

$$\begin{array}{rl}\Psi & \to T(k)\exp (-\mathrm{i}kx),x\to -\infty ,\\ \Psi & \to \exp (-\mathrm{i}kx)+R(k)\exp (\mathrm{i}kx),x\to +\infty ,\end{array}$$

$$S[u]:=\{R(k),-\infty <k<\infty ;p_{\alpha },q_{\alpha },\alpha =1,2,\dots ,A,\beta \in \mathbb{Z}\},$$

$$W(x,y)+K(x+y)+{\int }_x^{\infty }\mathrm{d}zW(x,z)K(z+y)=0,y>x,$$

$$K(x)=\frac{1}{2\pi }\int \mathrm{d}k\exp (\mathrm{i}kx)R(k)+\sum _{\alpha =1}^{\alpha }{q}_{\alpha }\Gamma \left[\frac{1}{\beta +1},p_{\alpha }(-x{)}^{\beta +1}\right],$$

$$\Gamma (a,z):={\int }_z^{\infty }\mathrm{d}s\exp (-s)s^{a-1},$$

$$S[u]=\{R(k)=0,-\infty <k<\infty ;p_1=p,q_1=q,A=1,\beta =1\},$$

$$W(x,y)+q\Gamma \left[\frac{1}{2},p(x+y{)}^2\right]+\zeta q{\int }_x^{\infty }\mathrm{d}zW(x,z)\Gamma \left[\frac{1}{2},p(z+y{)}^2\right]=0,y>x,$$

$$K_1(z+y):={\int }_z^{\infty }\mathrm{d}{z}^{\prime }K(z+z^{\prime })K(z^{\prime }+y).$$

$$-W(x,y)=K(x+y)-\zeta {\sigma }_{\mu }(x,y)-{\rho }_{\mu }(x,y),$$

$${\sigma }_{\mu }(x,y)$$

$${\rho }_{\mu }(x,y)$$

$$K_{\mu }(z+y)={\int }_z^{\infty }\mathrm{d}{z}^{\prime }{K}_{\mu -1}(z+z^{\prime })K(z^{\prime }+y).$$

$$\left|{\zeta }^{\mu -1}{\int }_x^{\infty }\mathrm{d}z{K}_{\mu }(z+y)K(x+z)\right|\le ϵ(|\zeta |ϵl{)}^{\mu -1}\Vert K(x+z)\Vert .$$

$$\begin{array}{rl}|{\sigma }_{\nu }(x,y)-{\sigma }_{\mu }(x,y)|& \le ϵ\left[\sum _{\kappa =\mu +1}^{\nu }(|\zeta |ϵl{)}^{\kappa -1}\right]\Vert K(x+y)\Vert \\ & \le ϵ\frac{(|\zeta |ϵl{)}^{\mu }}{1-|\zeta |ϵl}\Vert K(x+y)\Vert \\ & \le \delta ,\end{array}$$

$$|{\rho }_{\mu }(x,y)|\le ϵ(|\zeta |ϵl{)}^{\mu }\Vert W(x,z)\Vert ,$$

Full text

1. Introduction

The Schrödinger–Poisson system is relevant for several fields in physics and mathematics. There are cosmological models where the dark matter takes the form of a scalar field. Such field dynamics for galactic halos is described by a Schrödinger–Poisson system [1, 2, 3, 4, 5]. Hypothetical dense scalar field known as Boson stars are the relativistic counterpart of a galactic halo [6, 7, 8, 9, 10, 11, 12]. In the astrophysics realm, Schrödiger equation describes the dynamics of a scalar field while Poisson equation dictates the dynamics of the gravitational field. On the other hand, in quantum mechanics and semiconductor theories the Schrödinger–Poisson system models the interaction of charged particles in crystals [13, 14], e.g.. Schrödinger equation describes a particle dynamics under the presence of a electric field given by Poisson equation.

A convenient way to solve the Schrödinger–Poisson equations is by numerical computations [15]. Although there are several theoretical studies showing some of the general properties of the system, in particular its solitonic nature [16, 17, 18].

The first techniques to solve some celebrated solitonic equations were laid soon after Zabusky and Kruskal reported their observations on the Fermi–Pasta–Ulam puzzle [19] (even Seeger et al. noticed quite the same over a decade before [20]). Among these techniques, the spectral transform, originally introduced by Gardner et al. [21], is considered as the cornerstone of soliton theory given it is an extension of the Fourier transform to the nonlinear framework [22]. The theory became augmented by relevant enhancements and additional developments to the spectral transform [23, 24] as well as the introduction of alternative approaches founded on the Darboux and Bäcklund transformations [25, 26, 27, 28, 29]. Today, the length of these methods is such that these reach applications of indubitable interest in other areas of mathematical physics [30].

In this work we concentrate on the spectral transform in order to seek for exact, stationary solutions to the Schrödinger–Poisson equations. Foremost, our procedure is entirely heuristical and lies substantially on integral equations.

The paper is organized as follows: In section 2 approximated solutions are obtained through the homotopy analysis method [31]. This algorithm enable us to solve nonlinear integral equations, being an analogous path to the inverse scattering method where the linear integral equation of Gelfand–Levitan–Marchenko (GLM) raises (section 3). Finally, our results are summarized in section 4. The presented solutions are single–solitonic, i.e. are uniquely related to an individual bounded state.

1.1. Units and Notation

In order to simplify the writing, we use dimensionless quantities along the document, unless it is stated otherwise. In particular, we use geometrical units in the very first part of section 2 to introduce the Einstein–Klein–Gordon equations. From section 3, Schrödinger equation appears in its canonical form, where $\hbar=1$.

We try to avoid all functional dependences unless these are absolutely necessary because of ambiguity or lack of context. Also, the following notation is employed:

• $\mathrm{i}$ is the imaginary unit ($\sqrt{-1}$),
• $x,y$ are real independent variables,
• $\iota,\kappa,\mu,\nu$ are discrete (integer) indexes (quantities),
• $u=u(x)$ is an arbitrary test function,
• $u_{x}$ denotes the derivative of the function $u$ with respect to $x$,
• $R_{\mu\nu}$ is a tensor,
• $K(\cdot)$ is a kernel,
• $K_{\mu}(\cdot)$ is an iterated kernel,
• $\psi,\Psi$ are wave functions,
• $S[u]$ is the spectral transform of the function $u$,
• $\partial_{x}$ denotes a partial differential operator acting on the independent variable $x$,
• $\bigtriangleup$ stands for the Laplace operator,
• ${\bf H}$ is the Schrödinger operator,

unless a different thing is specified on the text.

2. First approximations

The Schrödinger–Poisson equations

[TABLE]

are the canonical representation of the Schrödinger–Newton equations originally introduced in [32], which, constitute the weak field limit of Einstein–Klein–Gordon equations

[TABLE]

here $R_{\mu\nu}$, $g$ and ${\bf T}_{\mu\nu}$ stand for the Ricci tensor, the metric and the energy–momentum operator, all in geometric units.

From now on we restrict our study to one spatial coordinate, namely $x$. Additionally, we will assume stationary solutions of the form $\psi(x,t)=u(x)\exp(\pm\mathrm{i}\omega t)$, where $\omega$ is a suitable real constant, $t$ represents the time and $u:\mathscr{X}\rightarrow\mathscr{Y}$ is, in general, a holomorfic function such that $\mathscr{X}$ and $\mathscr{Y}$ are open subsets in $\mathbb{C}$. Thence, the nonlinear system in (1) is reduced to

[TABLE]

having introduced a complex parameter $b$, for generality. The choice in the sign related to $\omega$ is physically irrelevant. Also, notice that $u$ and $\phi$ are time independent.

Broadly speaking, $u$ can be computed in the neighbourhood $|x-x_{0}|$ by Padé approximants of the form $[\mu/\mu+1]_{U}(x)$, with $U=u_{x}/u$ (furthermore details can be found in [33, 34]). Yet, first approximations of $u$ can be obtained by un-coupling (3).

Consider (3a), which differentiated two times with respect to $x$ reads

[TABLE]

by making the substitutions $\phi_{xx}\rightarrow|u|^{2}$, $(b\phi-\omega)u\rightarrow u_{xx}$ and solving for $\phi_{x}$, (4) becomes

[TABLE]

Yet again, consider (3a), by differentiating it with respect to $x$, solving for $\phi_{x}$ and equating with (5), we obtain an uncoupled scheme of the Schrödinger–Poisson equations:

[TABLE]

that possess just the same integral curves than (3).

In particular, if we consider $u$ as a real function to be regular for all real values of $x$ and to vanish asymptotically rapid, viz.

[TABLE]

only the first term in the right–hand side of eq. (6) does contribute significantly if compared to the others, consequently

[TABLE]

which general Urysohn form is

[TABLE]

here $C_{i}$ are constansts determined from boundary conditions and the function $F(x)$ equals to $c_{1}+c_{2}x+c_{3}x^{2}+c_{4}x^{3}+f(x)$, where $c_{i}$ stand for constants and $f=f(x)$ is a real valued function to be continuous in the sense of Lipschitz, to assure that equation (9) has a unique solution.

In order to approximate the solution of the nonlinear integral equation (9), hither we focus on the homotopy analysis method. That being said, assume that each individual term of the Urysohn form (9) has a kernel $K(x,y)\in\mathbb{R}$. Let $\Xi_{\mu}(u^{3})\in\mathbb{R}$ be the $\mu$–th order deformation constraint given by

[TABLE]

then, the homotopy–series solution to (9) is computed through

[TABLE]

subject to the initial guess $u_{0}=f$. The last expression is a pivotal trace in our quest for exact solutions. A way to obtain general exact solutions to (3) is by the inverse scattering method, which is the spanned version of (11). See our discussion in section 3. In our current analysis, while the kernel in (11) is separable, given the interval $x\in[\chi_{1},\chi_{2}]$, convergence of series (11) is assured iff

[TABLE]

is undeniably satisfied111Cf. [35, 36] for a comprehensive discussion.. Thus, at the hand of allowing $f=\mathrm{sech}(ax)$, with $a$ some real constant, it does imply that uniquely the $C_{2}$ and $C_{4}$ terms in (9) will behave as required by (7), so the constants $c_{1}=c_{2}=c_{3}=c_{4}=0$ and the solution to (9) is approximated in virtue of

[TABLE]

indeed, the programme (13) lead to the exact solution to (9), iff the series (11) does converge. Unfortunately, these integrals may become complicated, hence, in most cases, (13) should be enough for numerical computation of the first rough solutions.

Moreover, the fact that we have restricted our study to those solutions compelled to suffice (7), is only a feasible gimmick to work for solutions of accesible interpretation in the physics environment. Nonetheless, a different class might be considered for another purposes.

In the next section we study the generalization of the former blueprint in accordance with a parallel essay from soliton theory.

3. Devising exact solutions

Straightaway, looking at the stationary Schrödinger equation defined in all $x\in\mathbb{R}$

[TABLE]

the discrete part of the operator ${\bf H}=-\partial_{xx}+u$, consists of $A$ negative eigenvalues $k^{2}=-p^{2}_{\alpha}$, with $p_{\alpha}>0$, $\alpha=1,2,\ldots,A$. To each of these eigenvalues corresponds a unique solution $\Psi_{\alpha}$ impelled to satisfy the boundary condition:

[TABLE]

Whereas the continuum part, featured by all positive real values of the eigenvalue $k^{2}$, does typify the solution $\Psi$ to accomplish the asymptotic boundary conditions:

[TABLE]

where $T(k)$ and $R(k)$ stand for the transmission and reflection coefficients222To behold how both (15) and (16) are certainly consistent with (7)..

Having said all that, the extended spectral transform $S$ of the function $u$ is defined as

[TABLE]

where $q_{\alpha}$ are the normalization factors and $\beta$ is a discrete parameter related to the GLM equation

[TABLE]

through the complex valued function

[TABLE]

where

[TABLE]

is the so called incomplete gamma function.

We request the extended spectral transform of (3) as the specified set

[TABLE]

and, correspondingly, the GLM equation becomes

[TABLE]

with $\zeta\neq 0$ a parameter in pursuance of convenience.

Since the kernel in (21) is not separable, it is a matter of appositeness to solve (18) for $W(x,y)$ in terms of von Neumann series.

In that event, consider the initial definition,

[TABLE]

hence, after the GLM equation has been iterated $\mu>1$ times, there results

[TABLE]

where

[TABLE]

and the $\mu$–th iterated kernel

[TABLE]

As the latter prescribes, provided that the function $K_{\mu}$ suffices the Picard iteration, then it is assured the existence and uniqueness of the solution to (18) by Lindelöf’s apophthegm, thusly, the convergence of the von Neumann series (24) remains to be tested.

3.1. Analysis of convergence

The expansion (23) comprises that if $|K(x+y)|$ is bounded by a number $\epsilon\in\mathbb{R}$ in any closed interval in $x$ of lenght $l$, hence, inequity $|K(x+y)|\leq\epsilon$ implies that $K_{\mu}(x+y)$ is also bounded for all $\mu\geq 2$, that is, $|K_{\mu}(x+y)|\leq\epsilon(\epsilon\,l)^{\mu-1}$. Therefore, each term of (24a) satisfies

[TABLE]

In particular, when $|\zeta|<1$, the sequence $\{\sigma_{\mu}(x,y)\}$ of partial sums is a Cauchy sequence for any positive real number $\delta$ iff

[TABLE]

is held in the very limit $\mu\rightarrow\infty$. Quite in the same way for

[TABLE]

which means that $\rho_{\mu}(x,y)$ will vanish uniformly along the interval $l$ while $\mu\rightarrow\infty$. Thereupon, the sequence $\sigma_{\mu}(x,y)$ of continuous functions converges absolute and uniformly to the function

[TABLE]

given $\Xi(z,y;\zeta)$ as the dissolvent kernel, analogous to the higher order deformation constraint (10) in the ambiance of homotopy analysis studied in section 2.

Thus, the solution to (18) has the form

[TABLE]

and, on account of that, from (21) we obtain

[TABLE]

whilst the solution to equations (3) is computed through the formulae:

[TABLE]

In our particular case, we claim for a real valued, regular function $u(x)$ that vanishes asymptotically, exponentially,

[TABLE]

Inasmuch as this scenario holds, the reflection coefficient $R(k)$ is meromorphic in the Bargmann strip

[TABLE]

where there is a bijection between the poles of $R(k)$ and the discrete eigenvalues $k_{\alpha}=ip_{\alpha}$ [37, pp. 68–79]. For all these poles and all the discrete eigenvalues, the correspondence

[TABLE]

is prevailed.

On the other hand, as we have settled in (20), the reflection coefficient $R(k)$ is demanded to be zero for all $k$, thus, the last relationship is trivial and, whether or not our solutions are normalizable, it constitutes a relaxed constraint.

3.2. Sketching the solutions

In general, the dissolvent kernel $\Xi(z,y;\zeta)$ can be built as the quotient of two functions $P(z,y;\mu)$ and $Q(\mu)$ prescribed by the series expansion:

[TABLE]

with the proviso of initial data $\Lambda_{0}(z,y)=K(z+y)$ and $\lambda_{0}=1$. As for the coming terms $\mu\geq 1$, the function $\Lambda_{\mu}(z,y)$ is computed by means of the recursive relation:

[TABLE]

with

[TABLE]

recalling that $l\in\mathrm{Re}\,\mathscr{X}$.

After the algorithm (35)–(36) has been iterated one time, fixing $\zeta=1$ without loss of generality, the integral in (30) is explicitly calculated as

[TABLE]

this integral must converge to zero in the limit $\mu\rightarrow\infty$ of iterations since $\sigma(x,y)$ satisfies the Hölder condition for some positive real constant $\delta$ when becomes characterized as from (19). Then, as it decays monotonically with each iteration, the same reasons stated in 3.1 are concedable. As a result we drop out the second term from (30) and, therefore, the solution to (21) is simply

[TABLE]

Now, it is admitedly suggested that by the direct application of formulae (31) to (37), the solution to the system (3) has the general form

[TABLE]

where $E(x,x_{0})$ is a real valued, continuous function to be determined (as we will see in section 4), $x_{0}$ is a point in the neighbourhood of $x$ and $c_{1},c_{2}$ are integration constants.

Furthermore, the function $E(x,x_{0})$ must suffice (by itself) equation (14) and, whenever it vanishes faster than (32), namely

[TABLE]

the relationship (33) holds; consequently $u(x)$ to be normalizable is not a necessary condition and the factor $q$ is not relevant in our quest for exact solutions.

4. Detachable solutions

Let us consider the solution (38) and let us assume that $E_{1}=E_{1}(x,x_{0})$ and $E_{2}=E_{2}(x^{\prime},x^{\prime}_{0})$ satisfy by themselves the nonlinear system (3). We can look for separable solutions by starting with the stationary Schrödinger equation (14), from which the product $E_{1}E_{2}$ has to be solvable, i.e.

[TABLE]

such that both $\phi=\phi_{12}=\phi_{21}$ and $\phi=\phi_{1}+\phi_{2}=\phi_{2}+\phi_{1}$ must be fulfilled. The general solutions to the latter functional equation are $\phi=-a\ln|E|^{2}$ and $\phi_{\nu}=-a\ln|E_{\nu}|^{2}$, for some real constant $a$ and $\nu=1,2$. More than that, it is required from (38) that $E$ satisfies simultaneously the Poisson equation, that means

[TABLE]

In spirit of section 2, assume that $E$ is a continuous real function which behaves according to the law (7), hence (40) has a family of parametric solutions $E\rightarrow E(\cdot;\xi_{1},\xi_{2})$, namely,

[TABLE]

with $\xi_{1},\xi_{2}$ the real parameters which characterize the solutions. Notice that, in the very particular case $\xi_{1}=(2p)^{-1}\ln(2ap^{2})$ and $\xi_{2}=-(2p)^{-1}\ln(2ap^{2})$, the expression in (41) is reduced to

[TABLE]

Finally, with the aim of (38) and (41), and recalling the form of our initial ansatz $\psi(x,t)\rightarrow\psi(x,t;\xi_{1},\xi_{2})=u(x;\xi_{1},\xi_{2})\exp(\pm\mathrm{i}\omega t)$, we write down the general stationary family of solutions to the nonlinear system (3) as

[TABLE]

which are Schrödinger–Poisson, stationary, single, parametric solitons subject to the conservative dispersion relation:

[TABLE]

where we have neglected $\xi_{1},\xi_{2}$ and $x_{0}$ for a tight writing.

The wave function in (43) corresponds to a quiescent soliton centred at $x_{0}$, oscillating with speed $\omega$ and decaying according to the eigenvalue $p$ associated to a single bounded state.

In case $A\neq 1$ in (20), we would obtain the equivalent multiple–stationary soliton solutions once applied the procedure described above. This is not an easy task: the more bounded states, the more demanding to solve a system of $A$ GLM equations in a closed form. For this purpose, a numerical approach might be an alternative to achieve the results.

As a final remark, note that the number of discrete eigenvalues $A$ of the Schrödinger operator ${\bf H}$ in (14) is computed through a function $J=J(x)$ defined by the nonlinear equation

[TABLE]

subject to the boundary condition

[TABLE]

$c$ standing for any positive constant, then

[TABLE]

Certainly, $\omega$ is discretized since $a$ and $b$ are constants and $p$ is a discrete eigenvalue of ${\bf H}$, that is, because we have considered a single soliton ($A=1$) there is a single value for $\omega$. In consequence, for these particular case, (38) must satisfy (45) such that $1=\lfloor J(\infty)/\pi\rfloor$, just as it is. Properly, each multiple soliton solution to (1) has to fulfill (45) so the number of nonlinear superpositions (or Bianchi permutations) coincides with (47).

5. Discussion

We presented exact, stationary, parametric solutions to the Schrödinger–Poisson nonlinear system of partial differential equations. In the first part of our study, we showed how the homotophy analysis method suffices to integrate the system (1). However, we found that to compute the integrals in (13) is not straightforward, but instead, a limitation carried with this programme. As obvious, this approach might provide first numerical computations.

In the second part of our study, we found a family of exact solutions to the Schrödinger–Poisson system through the inverse scattering method. The inherent relation between these two techniques has been tangentially depicted.

In both cases, we addressed our discussion to the theory of integral equations, either nonlinear and linear. To solve these equations is not an easy task, nonetheless, the techniques come up with a deep understanding about the very nature of the problem in itself. At this point, we have not been able to seek for exact solutions by any algebraic method. Also, without counting the inverse scattering method, the mainstream techniques described in soliton theory are insufficient to solve our problem in a closed form.

Our solutions just cover the stationary configuration of the Schrödinger–Poisson system and are related to a single eigenvalue of the Schrödinger operator. To seek for closed solutions out of equilibrium is not a goal we have achieved yet. Quite the same for multiple Bianchi permutations, i.e. more than one eigenvalue (bounded state) of the Schrödinger operator. In that case, we should find something at the same level of a Bäcklund (or Darboux) transformation in order to induce a nonlinear superposition of bounded states.

Acknowledgments

The authors are glad to thank the valuable remarks and discussions by their colleagues at the Department of Physics at CINVESTAV, Mexico City.