An Approximation Method for the Scattering Data of One-Dimensional Soliton Equations under Arbitrary Rapidly-Decreasing Initial Pulses

TL;DR

This paper introduces a new approximation technique to predict the number of solitons generated by arbitrary initial pulses in one-dimensional soliton equations, avoiding complex integrations.

Contribution

The method provides a simple way to estimate soliton counts for various equations, demonstrated on the nonlinear Schrödinger equation with accurate results.

Findings

Accurately estimates soliton numbers without solving original equations

Shows good agreement with numerical simulations

Applicable to a wide range of one-dimensional soliton equations

Abstract

We present a novel approximation method that can predict the number of solitons asymptotically appearing under arbitrary rapidly decreasing initial wave packets. The number of solitons can be estimated without integration of the original soliton equations. As an example, we take the one-dimensional nonlinear Schrodinger equation and estimate the behaviors of the scattering amplitude in detail. The results show good agreement compared with those obtained by direct numerical integration. The presented method is applicable to a wide class of one-dimensional soliton equations.

Equations

$${\Psi }_x$$

$${\Psi }_t$$

$$i{\psi }_t=-{\psi }_{xx}-2|\psi {|}^2\psi .$$

$$S$$

$$T$$

$$\psi(x,0)\to 0,\quad\mbox{as $|x|\to\infty$.}$$

$$\varphi (x;\xi )\to \left(\begin{array}{c}{e}^{-i\xi x}\\ 0\end{array}\right),\bar{\varphi }(x;\xi )\to \left(\begin{array}{c}0\\ e^{i\xi x}\end{array}\right),$$

$$\chi (x;\xi )\to \left(\begin{array}{c}0\\ e^{i\xi x}\end{array}\right),\bar{\chi }(x;\xi )\to \left(\begin{array}{c}{e}^{-i\xi x}\\ 0\end{array}\right),$$

$$\begin{array}{r}\varphi (x;\xi )=a(\xi )\bar{\chi }(x;\xi )+b(\xi )\chi (x;\xi ),\\ \bar{\varphi }(x;\xi )=\bar{a}(\xi )\chi (x;\xi )-\bar{b}(\xi )\bar{\chi }(x;\xi ).\end{array}$$

$$\phi(x;\xi)=\begin{pmatrix}a(\xi)e^{-i\xi x}\\ b(\xi)e^{i\xi x}\\ \end{pmatrix}\quad\mbox{as $x\to+\infty$.}.$$

$$K_1^{\ast }(x,y;t)-F(x+y;t)-{\int }_x^{\infty }{K}_2(x,z;t)F(z+y;t)dz=0,$$

$$K_2^{\ast }(x,y;t)+{\int }_x^{\infty }{K}_1(x,z;t)F(z+y;t)dz=0,$$

$$F(x;t)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\frac{b(\xi )}{a(\xi )}{e}^{i\xi x-4i{\xi }^{2}t}d\xi -i\sum _{l=1}^M\frac{b({\zeta }_l)}{a^{\prime }({\zeta }_l)}{e}^{i{\zeta }_{l}x-4i{\zeta }_{l}t}.$$

$$\psi (x,t)=-2i{K}_1^{\ast }(x,x;t).$$

$${\Psi }_x=S\Psi ,S=\left(\begin{array}{cc}-i\xi & i{\psi }^{\ast }\\ i\psi & i\xi \end{array}\right),$$

$$I_j:x_j\le x<x_{j+1}(j=1,\dots ,N),$$

$${\psi }_j(x)=\{\begin{array}{ll}{V}_j& x\in I_j,\\ 0& x\notin I_j.\end{array}$$

$$\psi (x,0)$$

$$=\{\begin{array}{ll}{V}_j& (x\in I_j,j=1,2,\dots ,N),\\ 0& \text{mbox}(otherwise).\end{array}$$

$${\Psi }_x=S_j\Psi ,S_j=\left(\begin{array}{cc}-i\xi & i{V}_j^{\ast }\\ i{V}_j& i\xi \end{array}\right).$$

$$\Psi (x,\xi )=\left(\begin{array}{c}{\Psi }_1(x,\xi )\\ {\Psi }_2(x,\xi )\end{array}\right),{\Psi }_0(x,\xi )=\left(\begin{array}{c}{e}^{-i\xi x}\\ 0\end{array}\right).$$

$$g(x-x^{\prime })=\left(\begin{array}{cc}{e}^{-i\xi (x-x^{\prime })}\theta (x-x^{\prime })& 0\\ 0& e^{i\xi (x-x^{\prime })}\theta (x-x^{\prime })\end{array}\right),$$

$$U(x)=\left(\begin{array}{cc}0& i{\psi }^{\ast }(x,0)\\ i\psi (x,0)& 0\end{array}\right).$$

$$\Psi (x,\xi )={\Psi }_0(x,\xi )+{\int }_{-\infty }^{\infty }d{x}^{\prime }g(x-x^{\prime })U(x^{\prime })\Psi (x^{\prime },\xi ).$$

$$\begin{array}{rl}{\Psi }_1(x,\xi )=e^{-i\xi x}+& e^{-i\xi x}{\int }_{-\infty }^{\infty }d{x}^{\prime }d{x}^{\prime \prime }\theta (x-x^{\prime })e^{2i\xi x^{\prime }}i{\psi }^{\ast }(x^{\prime },0)\\ & \cdot \theta (x^{\prime }-x^{\prime \prime })e^{-2i\xi x^{\prime \prime }}i\psi (x^{\prime \prime },0)+(\mathrm{even}).\end{array}$$

$$a(\xi )e^{-i\xi x}=\underset{x\to \infty }{\lim }{\Psi }_1(x,\xi )=\underset{x\to \infty }{\lim }{\varphi }_1(x,\xi ).$$

$$v(x)=i{e}^{2i\xi x}{\psi }^{\ast }(x,0),$$

$$v^{\ast }(x)=i{e}^{-2i\xi x}\psi (x,0).$$

$$a(\xi )=1+\underset{x\to \infty }{\lim }A(x;\xi ),$$

$$\begin{array}{rl}A(x;\xi )={\int }_{-\infty }^{\infty }d{x}^{\prime }{\int }_{-\infty }^{\infty }d{x}^{\prime \prime }& \theta (x-x^{\prime })v(x^{\prime })\\ & \cdot \theta (x^{\prime }-x^{\prime \prime })v^{\ast }(x^{\prime \prime })+(\mathrm{even}).\end{array}$$

$$a(\xi )=1+\underset{n\to \infty }{\lim }A(x_n;\xi ).$$

$$A(x_j;\xi )=\sum _{k,l=1}^N\theta (x_j-x_k)v(x_k)\theta (x_k-x_l)v^{\ast }(x_l)(\Delta x{)}^2+(\text{even}),$$

$$v(x_k)=i{e}^{2i\xi x_k}{V}_k^{\ast },v^{\ast }(x_l)=i{e}^{-2i\xi x_k}{V}_l.$$

Full text

1 Introduction

Soliton theory is one of the key concepts in mathematical physics. It is particularly important for continuous systems owing to the fact that a number of intriguing phenomena across various fields of physics can be systematically described in this common language of solitons. The fields where solitons play prominent roles include high energy physics[1, 2, 3, 4, 5] and gravitational physics[6]. As examples of physics in more easily accessible energy scales, we can cite continuum mechanics of fluids or plasmas[7, 8, 9, 10, 11, 12] and nonlinear optics[13, 14, 15, 16, 17]. In the context of condensed matter or low temperature physics, the macroscopic dynamics of Bose–Einstein condensed (BEC) systems is known to obey a kind of nonlinear Schr̈odinger equation (NLSE) with external potentials[18, 19, 20, 21]. These rich applications means that this mature subject is still an active area of research.

Soliton equations are nonlinear partial differential equations often discussed in the context of integrable systems. Although some subtleties exist in rigorously defining integrability for systems with infinite degrees of freedom, it is sufficient for our purpose to point out the applicability of the inverse scattering transformation (IST) method. This prominent property is shared by all soliton equations[22, 23, 24, 25, 26, 27].

By virtue of the IST method, soliton equations are decomposed into a pair of linear simultaneous differential equations in the following manner:

[TABLE]

where the quantities $S$ and $T$ are matrices or operators including the unknown functions of the original equation and an additional parameter called the spectral parameter. The wave function $\Psi$ represents an auxiliary field satisfying appropriate boundary conditions. Equation (1a) takes the form of an eigenvalue problem or a stable scattering problem where the initial condition plays the role of the scattering potential. Equation (1b) governs the time evolution of the scattering data whose initial values are determined by Eq. (1a).

Utilizing the obtained scattering data, we can construct exact solutions to the original equations at an arbitrary time. It is known that the original initial value problems can be solved under a fairly wide class of initial conditions. A class of rapidly decreasing functions is one of the most important classes of initial conditions for many soliton equations. Needless to say, it is almost impossible to prepare ideal pure solitons in actual experiments and hence it is more realistic to consider more general initial conditions rather than pure solitons. For many soliton equations, initial value functions should be rapidly decreasing, and this should be the only requirement that we can reasonably impose. One of the important quantities characterizing the initial conditions is the number of remaining solitons at $t\to\infty$. Since our target time domain lies in the asymptotic future and because of the diffusing nature of wave packets, the need for subtle treatment of boundary conditions arises. Therefore, a theoretical framework not depending on direct numerical integration will be useful and can serve to endorse the validity of numerical calculations.

To this end, we introduce a novel approximation method that can predict the number of solitons in the final state without direct numerical integration. The purpose of this paper is to present this method with some illustrative applications. This method provides a way to approximately obtain the scattering data for Eq. (1a) by discretizing the space coordinate and it is different from the method we used in our previous work[28]. The method presented in this paper is intended to overcome the drawbacks of our previously adopted semi–analytic method.

In contrast to the above treatments, the new approximation method presented in this paper is purely numerical and can be applied to realistic smooth functions for which sufficient fineness of the discretization pitch is required.

This paper is organized as follows. In sect.2, we present mathematical preliminaries by taking the NLSE as an example. Section 3 is the central part of this paper and we develop the method for detailed analysis of the initial value problem. In sec.4, we show some illustrative applications and present a brief model of a collision experiment with attractively interacting BEC solitons[29, 30, 31]. We compare the results obtained by the presented method with those from direct numerical integration to check the consistency between them. The final section is devoted to a discussions and concluding remarks.

2 Mathematical Preliminaries

2.1 Inverse scattering transformation

To be specific, we consider the NLSE with a self–focusing interaction throughout this paper,

[TABLE]

We assume that the initial wave packet is a rapidly decreasing function of $x$ hereafter. The auxiliary linear equations (1) for the NLSE are derived by introducing the two–element vector function $\Psi$ and 2$\times$2 matrices $S$ and $T$ as

[TABLE]

where $\xi$ is the time–independent spectral parameter. Many soliton equations can be formulated using the 2$\times$2 matrices in the framework of the IST method and they are classified by the way the spectral parameter $\xi$ appears in these matrices. The NLSE belongs to the Ablowitz–Kaup–Newell–Segur (AKNS) system[32], and other soliton equations belonging to the AKNS system can be discussed in a similar way.

An important step of the IST method is to analyze the spatial Eq. (1a) as a stable scattering problem, whose potential term is given by the initial condition $\psi(x,0)$. The wave function $\Psi$ and the spectral parameter $\xi$ correspond to the eigenfunction and eigenvalue, respectively. We assume the rapidly decreasing boundary condition for $\psi(x,0)$ as already mentioned:

[TABLE]

Equations (1a) and (3a) under the above boundary condition completely define the Zakharov–Shabat (ZS) problem[33] for the self–focusing NLS equation. With this boundary condition, each element of the wave function $\Psi$ must become a plane wave in the limit of $x\to\pm\infty$. As the fundamental solutions, we can select two sets of functions ${\{\phi,\bar{\phi}\}}$ and ${\{\chi,\bar{\chi}\}}$ called the Jost functions, which satisfy the boundary conditions

[TABLE]

The Jost functions are related to each other as

[TABLE]

The coefficients are members of the scattering data, and $a(\xi)$ can be analytically continued to the upper half plane $\mathop{\rm Im}\xi>0$.

From Eqs. (5) and (6), we can see that the Jost function $\phi(x;\xi)$ satisfies the asymptotic form

[TABLE]

Once the ZS problem is solved, the time–evolved wave function is easily obtained by Eq. (1b), and the solution under the initial value is formally constructed by virtue of the Gel’fand–Levitan–Malchenko (GLM) equation[34]:

[TABLE]

[TABLE]

where the scattering data appear through the function

[TABLE]

The solution of the NLSE is given by the kernel function $K_{1}$ as

[TABLE]

In the above equation, $\xi$ stands for the continuous real eigenvalue and $\eta$ are the discrete eigenvalues, which are the zeros of $a(\xi)$. The GLM equation clearly states that the number of discrete eigenvalues determines the number of asymptotically generated solitons. When the function $a(\xi)$ has $M$ simple zeros $\xi=\xi_{1},\xi_{2},\ldots,\xi_{M}$ on the upper half plane of $\xi$, $M$ solitons appear asymptotically and each $\xi$ determines the characteristics of each soliton.

2.2 Discretization of the initial wave packet

As explained in sect.2.1, we need to know $a(\xi)$ to extract the desired information of the final state. We find that this is equivalent to calculating $\phi(x;\xi)$ at $x\to+\infty$ under the boundary condition Eq. (5a).

In the ZS problem for the NLSE

[TABLE]

the function $\psi$ is the initial value of the unknown function $\psi(x,t)$. The ZS problem for the NLSE under the initial condition $\psi(x,0)=A\mathrm{sech}(x)$, where $A$ is a constant, has been extensively studied[35]. The explicit forms of the wave function $\Psi$ are written by the hypergeometric function and the number of remaining solitons is given by the maximum integer $n$ that satisfies $n\leq A$. Nevertheless, the ZS problems for more general initial conditions are difficult to solve, and it is hardly possible to predict how many solitons remain in the final state.

The major difficulty in analyzing Eq. (12) for general initial conditions originates from the fact that $\psi(x,0)$ depends on the coordinate $x$. In order to deal with this difficulty, let us split the support of $\psi(x,0)$ into many small intervals:

[TABLE]

and approximate $\psi(x,0)$ so that it takes a constant value in each interval. We introduce a set of functions $\psi_{j}$:

[TABLE]

The initial value $\psi(x,0)$ is now approximated as

[TABLE]

Assuming $\psi(x,0)$ belongs to the class of rapidly decreasing functions, we can approximately consider that $\psi(x,0)$ has a compact support. Within each interval, Eq. (12) reads as

[TABLE]

3 New Approximation Method

Now, let us develop our method with the discretization scheme. First, we introduce a two–component vector function $\Psi$:

[TABLE]

The 2$\times$2 matrix Green function $g(x-x^{\prime})$ is defined using the step function $\theta(x)$ as

[TABLE]

and the scattering potential in the matrix form is denoted by

[TABLE]

The starting point is to reformulate the ZS eigenvalue problem in Eq. (12) as Fredholm integral equations of the second kind:

[TABLE]

Needless to say, the solution $\Psi(x,\xi)$ asymptotically coincides with the Jost function in Eq. (7) at the limit $x\to+\infty$. By iteration, we can construct the perturbative solution to Eq. (21). By considering the first element of the two–component vector $\Psi(x,\xi)$, we can write down following Neumann series:

[TABLE]

According to Eq. (7), the transmission coefficient $a(\xi)$ can be obtained by taking large–$x$ limit:

[TABLE]

Note that only terms of even orders contribute to the transmission coefficient $a(\xi)$. This is because the number of iterative operations is equal to the number of scattering events that invert the direction of the incident plane wave $\Psi_{0}(x,\xi)$.

Similarly, by considering the second element of the two–component vector $\Psi(x,\xi)$, we can obtain the reflection coefficient $b(\xi)$, to which only terms of odd order contribute.

For later use, it is convenient to redefine the potential term to include the phase factor $e^{2i\xi x}$,

[TABLE]

[TABLE]

By using the above notations, the expression for the transmission coefficient $a(\xi)$ can be rewritten as

[TABLE]

We now consider the replacement of the integral in Eq.  (26) by a matrix operation of finite size. Although the reduced expression for $a(\xi)$ becomes an approximate one, it is more suitable for numerical calculation. We assume that the initial wave packet is distributed in a finite range of width $L$. In order to carry out the integral and the limit in Eq. (26), we use the discrete variables $x_{j}$ ($j=1,2,\ldots,N+1$) defined in Eq. (13), and we have the following approximate expression for $a(\xi)$:

[TABLE]

By using these discrete variables, the value of $A(x;\xi)$ at discrete values of $x$ can be expressed by a finite sum as

[TABLE]

where $\Delta x$ is the spacing of the discrete variables, given as $\Delta x=L/N$. Since the quantites $v(x_{k})$ and $v^{*}(x_{l})$ are defined in the intervals in Eq. (13), they are expressed by $V_{j}$ in Eq. (15) as

[TABLE]

We introduce matrices $V$ and $V^{*}$ as

[TABLE]

We also define a step function matrix $G$ whose $(i,j)$–element is given by $\theta(x_{i}-x_{j})\Delta x$. The matrix $G$ is a lower–triangular matrix, and its explicit form is

[TABLE]

The size of the initial wave packet should be sufficiently large so that the rapidly decreasing tails of the initial wave packet become negligible. To meet the convergence condition for the Neumann series, we can choose a sufficiently large number of small intervals so that $|V_{i}|L/N$ is always smaller than unity for every index $i$.

From the matrices $V$, $V^{*}$, and $G$, we find that

[TABLE]

This means that $A(x_{j};\xi)$ is given by the sum of the elements in the $j$th row of the matrix

[TABLE]

If we introduce the column vector \mbox{\boldmath{e}}=(1,\ldots,1)^{\rm T} and write the $j$th row of $W$ as \mbox{\boldmath{w}}_{j}, we obtain

[TABLE]

Since the scattering coefficient is given by Eq. (27) and the initial pulse is located in the region $x\leq x_{N+1}$, we can derive an approximate expression for $a(\xi)$ as

[TABLE]

This is equivalent to the expression

[TABLE]

Many useful algorithms for obtaining the inverse matrix are available and we can calculate Eq. (36) easily.

Theoretically, since the fact the matrix product of $G$ and $V$ is still a lower–triangular matrix, we can explicitly write down the formula for $a(\xi)$ as

[TABLE]

where the following notations are used:

[TABLE]

Although the factor $1/N^{2}$ exists in Eq. (38), the summation takes values between the first order and second order of the quantity $1/N$ as a whole.

4 Applications

We show some specific applications of the presented method. First, we calculate the transmission coefficient $a(\xi)$ for two different types of single–peak initial conditions, both of which are rapidly decreasing smooth functions. Through this simple example, the validity and effectiveness of our method is corroborated. Secondly, we consider two superposed wave packets whose centers are sufficiently separated so that the resultant initial packet forms a double–peaked profile. The phases of the two constituent pulses are varied to investigate the effect of the relative phases on the asymptotically generated states. The prediction obtained from our time–independent method is compared with the result from direct numerical integration of the original NLSE.

4.1 Sech–type and Gauss–type wave packets

We consider two typical initial conditions for the NLSE in Eq. (2), one of which is taken to be the pure one–soliton

[TABLE]

Both functions are shown in Fig. 1. Their profiles are similar. They also share the same $L^{2}$–norm and have almost the same full width at half maximum (FWHM). Within the FWHM, the difference between these two functions is no more than 4% of their maximum value of 2 at the origin. We take the values of $L$ and $N$ of $-3.0\leq L\leq 3.0$ and $N=1200$ and all the tails of these functions are truncated at $x=\pm 3.0$. We calculate the transmission coefficient $a(\xi)$ for these truncated wave packets using the method developed in the previous section.

The values of $|a(\xi)|$ for real $\xi$ under the initial conditions in Eqs. (39a) and (39b) are shown in Fig. 2. As expected from the pure spliton condition, the transmission coefficient $a(\xi)$ for Eq. (39a) is approximately unity and flat over a wide range of the spectral parameter, $0.0\leq\xi\leq 5.0$. By contrast, the plot of the transmission coefficient $a(\xi)$ for Eq. (39b) is not a constant. This indicates that the Gauss–type wave packet can never become a pure soliton. In this simple example, our method successfully shows the well–known result, which has been given by analytic proof. In Fig. 2, $|a(\xi)|$ is smaller than unity. This is partially because of the truncating operation that we applied. What matters more is the insufficient value of $N$. We next confirmed the approach of $|a(\xi)|$ to unity when $N$ was increased. In Fig. 3, we show plots of $|a(\xi)|$ for various $N$, where all the plots were calculated for the same initial wave packet Eq. (39a). In Figs. 4 and 5, visualized two–dimensional (2–D) plots of the matrix function $\mathrm{log}|1/(1-GVGV^{\ast})|$ for Eqs. (39a) and (39b) under $\xi=0.0$ are respectively shown. These matrices are both lower–triangular. The last rows of them, which are used for the summing up, are simultaneously shown in Fig. 6. The sum of the last $N$–th row for Eq. (39a) appears to take larger values around a narrow neighborhood of the origin $x=0.0$. For Eq. (39b), the distribution is more widely spread. This reflects the feature of the pure soliton, that is, certain information is maximally concentrated. This kind of numerical experiment brings some insight that enables us to characterize the geometrical profile of pulse like functions.

4.2 Double–peaked pulses

Next, we consider some examples of double–peaked initial wave packets. This type of initial condition is given by a superposition of two wave packets, each of which has a single peak that is appropriately apart from the other. There can be relative phases $\Delta\phi$ between these constituent wave packets. We show three typical examples of simulations for different values of $\Delta\phi$. Throughout this subsection, we set $N=600$ and the increment of $\xi$ to be 0.01 along the real and imaginary axes.

4.2.1 Case I: $\Delta\phi=0$

We consider the initial wave packet

[TABLE]

No relative phase exists between the two constituent pulses in this case. The profile of this wave packet is shown in Fig. 7. We can see that it has two peaks at $x=\pm 1.0$. The absolute value $|a(\xi)|$ on the upper half plane of $\xi$ calculated by the method presented in the previous section is plotted in Fig. 8. The same quantity is also depicted three–dimensionally (3–D) in Fig. 9. The distribution is symmetric about the imaginary axis as expected, and a zero exists around $\xi=0.31i$. This means that only one soliton remains asymptotically. To verify this prediction, we numerically integrated the original NLSE. The absolute squares of the solution $\psi(x,t)$ at $t=5.0$, $15.0$, and $25.0$ are shown in Fig. 10, and we confirmed that only one pulse was generated around the origin.

4.2.2 Case II: $\Delta\phi=\pi$

We set the initial wave packet to be

[TABLE]

This time, the relative phase $\Delta\phi$ is taken to be $\pi$. The profile of this wave packet is shown in Fig. 11. We can see that it has two peaks with opposite convexity at $x=\pm 1.0$. The absolute value of $a(\xi)$ for this case is plotted two–dimensionally in Fig. 12 and three–dimensionally in Fig. 13. On the upper half plane of $\xi$, these plots are again symmetric about the imaginary axis, and two zeros appear around $\xi=\pm 0.62+0.17i$. This means that two solitons will remain in the final state. We confirmed this prediction by direct numerical integration. As shown in Fig. 14, we can observe two pulses symmetrically located around the origin as expected. The magnitudes of the two solitons are smaller than that of the single soliton generated in Case I. This is because the two pulses in Fig. 14 are wider than the single pulse shown in Fig. 10. Actually, the heights of these two pulses remained at the level shown in Fig. 14 even if we extended the termination time of the numerical integration.

4.2.3 Case III: $\Delta\phi=\frac{\pi}{2}$

In this case, we take initial wave packet as

[TABLE]

The relative phase $\Delta\phi$ is set to $\pi/2$, which corresponds to the situation where one pulse approaches the other, then the two pulses collide. A 2–D plot of $|a(\xi)|$ is shown in Fig. 15 and a 3–D plot of $|a(\xi)|$ is shown in Fig. 16. For this case, both plots are no longer symmetric. We can find two zeros of $a(\xi)$, one zero around $\xi=0.33+0.27i$ and the other around $\xi=0.93+0.01i$. For this initial condition, we also executed a numerical integration and observed two major pulses, a larger one and a smaller one, as shown in Fig. 17. Their centers are located asymmetrically about the origin.

In the above examples, the number of finally remaining solitons is at most two, but three or more solitons are generated if we choose suitable initial conditions. In fact, we showed that two box–type pulses with suitable height and separation develop to generate three solitons[28]. Concerning the number of solitons in the final state, the presented method has no limitations in principle, as long as $N$ is taken to be sufficiently large. We have confirmed that our previous result is also reproduced by the method presented in this paper.

An actual experiment on soliton collision has been performed, where attractively interacting BECs were exploited[29]. In the experiment, a confinement trap in the x–direction was applied and was not completely removed during the collision process. Owing to the existence of this trap, the geometry was not completely flat. The trapping potential is usually well–approximated by a quadratic function. It is known that the AKNS–ZS formulation for the NLSE can be extended to a situation where such a potential term exists[36].

As long as we choose suitable initial conditions, by transforming a spatial variable x to a new one and reconstructing a new scattering potential from a given initial wave packet, we can use the same procedure as that demonstrated in previous sections. Nevertheless, it is unclear whether all initial conditions with a rapidly decreasing property work well in this extended scheme. There is a possibility that this method imposes additional restrictions on the class of initial conditions, which leads to a loss of generality. To avoid this complication, we neglected the effect of the trapping potentials in our illustration.

In the BEC experiment referred to above[29], the effects of the relative phase on the collision process were examined and it was reported that the behaviors of the two condensate pulses were different in the in–phase and out–of–phase cases. To extract and discuss the role of the relative phase between condensates, our simplified models seem to be suitable and sufficient. In this sense, it can be reasonabley stated that these examples correspond to real physics and have physical interest.

5 Summary and Discussion

In this paper, we have presented a novel approximation method for obtaining the scattering data for the ZS problem reduced from the one–dimensional NLSE under general rapidly decreasing initial conditions. Using the obtained scattering data, the numbers of asymptotically generated solitons have been estimated for initial conditions whose profiles are not limited to pure solitons. The derived results are consistent with those of direct numerical integration. These results include brief models for a collision experiment on BEC solitons. A similar experiment on soliton collision in a homogeneous space is expected to be more easily conducted in the field of nonlinear optics. Because of its ubiquitous appearance and importance in physics, we have limited our illustration to the NLSE. The method developed in this paper is essentially applicable to a wide class of soliton equations that can be formulated as the AKNS, Kaup–Newell,[37] and Wadati–Konno–Ichikawa[38] systems. Other soliton equations of importance that can be fit into this format include the Korteweg–de Vries (KdV), modified KdV, sine–Gordon, and derivative NLS equations.

Since the NLSE is a Hamilton system, an action variable $J(\xi)$ and an angular variable $\theta(\xi)$ exist, with which the Hamiltonian can be written. Excluding the zeros, they are expressed by the scattering data as:

[TABLE]

[TABLE]

Using the technique developed in this paper, we have obtained a new numerial method of investigating the geometric surfaces defined by these action–angle variables. Although we have shown only the absolute values $|a(\xi)|$ for visualization, the original complex geometry of the scattering data in the upper half plane of $\xi$ is expected to include richer information. Our ultimate goal is to find a way to extract physical information through the investigation of these geometrical entities. It is desirable that we can define some characteristic quantities concerning the geometry so that they correspond to actual physical information. This task is left as future work.

{acknowledgment}

The authors express their sincere gratitude to Professor Ralph Willox of the University of Tokyo and Professor Ken–ichi Maruno of Waseda University for their interest in this work and stimulating discussions. One of the authors (H. F.) thanks Utsunomiya University for offering opportunities for fruitful discussions.