MatterWave Fields for DoubleSlit Atom Interferometry: Variational Versus Exact Solitons
\abstA major challenge in the theoretical modeling of doubleslit interferometry involving matterwave fields is the appropriate waveform to be assigned to this field. While all the studies carried out to date on this issue deal with variational fields, experiments suggest that the optical field is generated by splitting a singlehump Bose–Einstein condensate into two spatially and temporally entangled pulses indicating the possibility of fully controlling the subsequent motion of the two output pulses. To probe the consistency of variational and exact soliton solutions to the field equation, we solve the Gross–Pitaevskii equation with an optical potential barrier assumed to act as a beam splitter, while including gravity. The exact solution is compared with the two most common variational wavefunctions, namely, the Hermite–Gaussian and supersech modes. From numerical simulations, evidence is given of the exact solution as being the most appropriate matterwave structure that provides a coherent description of the generation and spatiotemporal evolution of matterwave optical fields in a hypothetical implementation of doubleslit atom interferometry.
1 Introduction
Boundsoliton pairs provide interesting transmission channels for multiplexed highintensity pulse trains [1, 2, 3, 4, 5, 6]. Among them, bisolitons have attracted a great deal of attention following their prediction [7, 8, 9, 10, 11, 12, 13, 14] and observation [15, 16] in several distinct optical media. In optical fibers, for instance, such structures originate from splitting [17, 18] a singlepulse optical soliton using a quadratic chirp, which leads to a soliton molecule whose intensity profile exhibits two temporal and spectral peaks with a finite phase difference between them [8, 16]. Similar objects have been predicted and observed in Bose–Einstein condensate (BEC) systems [19, 20, 21], where they result from a beam splitter acting on a bright matterwave pulse, which creates two copropagating matterwave pulses.
Bisolitons are of fundamental importance in BECs since they provide appropriate wave structures needed for the experimental achievement of an atomic equivalent of doubleslit holography [22, 23]. In such an experiment, we can imagine two spatially and temporally entangled matterwave pulses generated from the splitting of a single condensate into two atomic populations, namely in nearly degenerate hyperfine levels (see for example ref. [24]). Generally a faroff resonant laser barrier provides an ideal scheme for such splitting [22, 23]. In theory, this most often has been described by considering a single matterwave pulse passing through a potential barrier created by a laser field, with only a fraction of the macroscopic atomic system allowed to cross the barrier while the other fraction remains trapped on the other side of the barrier.
A pioneering model for such process was proposed in ref. [25], in terms of a single macroscopic atomic population moving in a doublewell optical field of a finite barrier. In this model, the system dynamics leads to two unbalanced population fractions from a single initial atomic population trapped in the two wells surrounding the potential barrier, with one fraction in the left well and the other fraction in the right well, after tunneling through the finite barrier. In the approach of ref. [25], however, the focus was on the quantum tunneling dynamics of atoms, particularly the associated Josephson effect related to the phase difference between the wavefunctions of the two condensates. Nevertheless, having at hand an analytical expression for the full condensate wavefunction can also be useful for experiments. This latter problem can be addressed by considering the Gross–Pitaevskii (GP) equation with an antiharmonic (i.e., a repulsive harmonic) optical potential representing the laser barrier, which is actually equivalent to the nonlinear Schrödinger equation (NLSE) describing the propagation of pulse signals in an optical fiber with a quadratic chirp [8, 9, 10, 12, 13, 14, 26, 27].
In recent works dealing with matterwave solitons within the framework of the GP equation with a harmonic optical potential, emphasis has mostly been on variational considerations, where the macroscopic condensate wavefunction is approximated by the Hermite–Gaussian mode or a supersech pulse [28, 29, 30]. Still, the variational treatment is based on perturbation theory, both in the choice of the trial solution (despite its localized pulse shape, the Hermite–Gaussian mode is not a solution to the GP equation) and in the formulation of the time evolution of the variational soliton’s collective coordinates [31].
In this work, we address the issue by considering two distinct variational solutions, namely, the supersech and Hermite–Gaussian modes [19, 28, 29, 30] on one hand, and the exact onesoliton solution to the associated GP equation obtained by a nonisospectral inversescattering transform (NIST) method [32] on the other hand. A gravitational potential is included with the aim of taking into consideration the possible relative acceleration of one pulse in the postcreated bisoliton with respect to the other due to the freefall motion of atoms in the condensate fraction exposed to gravity [33, 34, 35, 36]. In effect, because of the finite mass of atoms, gravity turns out to be an additional fundamental ingredient in matterwave interferometry as previously emphasized (see for example refs. [22, 23]).
2 Exact OneSoliton Solution
Consider a system of weakly interacting identical atoms of mass . Because of their interactions, the atoms form an ensemble falling freely under gravity but also experiencing a repulsive force due to a potential barrier of finite height. In the meanfield picture, the quantum states of such a macroscopic system of bosons can be represented by a wavefunction , where represents the direction of the free fall of atoms and is the time variable. The GP equation governing the spatial and temporal evolution of is given by
(1) 
where is the mean–field interatomic interaction coefficient related to the swave scattering length. is the total external potential for the nonlinear Schrödingertype Eq. (1), it is given as the sum of the repulsive harmonic potential and the linear potential accounting for the gravitational field
(2) 
with the acceleration of gravity and a characteristic frequency. To obtain the exact solution to Eq. (1) using the NIST method, it is useful to introduce the following dimensionless variables
(3)  
(4)  
(5) 
With these new variables, the GP equation in Eq.(1) becomes
(6) 
which is simply the selffocusing NLSE with a nonlocal (i.e., spatially varying) external potential . When and , Eq. (6) describes the macroscopic dynamics of a freefalling BEC. The same equation was proposed to model the dynamics of plasmas in a gravitational field, and its exact one and nsoliton solutions have been obtained [37] by means of the InverseScattering Transform (IST). On the other hand, when but , Eq. (6) becomes the NLSE with an external spacevarying repulsive harmonic potential. Exact soliton solutions for this latter physical context have also been obtained, using an improved IST technique that overcomes the nonisospectral character of the associated Laxpairtype eigenvalue problem [38].
When both and are nonzero, the perturbed NLSE Eq. (6) gives rise to an interesting problem in mathematical physics whose solution requires a subtle combination of the methods proposed independently for the two particular cases: but [38, 39, 40, 41], and with [37]. In the general picture of the IST [42, 43] the first step is the construction of two pairs of coupled linear eigenvalue equations associated with the perturbed NLSE Eq. (6). Thus, let us consider the two following coupled pairs of eigenvalue problems [32]
(7) 
(8) 
where and are spacetimedependent functions that form a twocomponent spinor with the associate eigenvalue, while
, and are unknown functions that can be computed from the condition of compatibility of the coupled pairs of equations.
The quantity , which is precisely the field function of our main interest, is here regarded as a scattering potential for the Laxpair type
eigenvalue problem. In the spirit of the IST [42, 43], exact solutions to the eigenvalue problems in Eqs. (7) and (8) follow from the scattering matrix for a specific choice of the scattering potential.
For a localized (i.e., pulseshaped) solution, we will require that the scattering potential vanishes asymptotically as at any propagation time. Such
asymptotic behavior is only possible for a judicious choice [44] of the initial
scattering potential , and from this choice the time evolution of the scattering data can easily be computed
via GelfandLevitanMarchenko [42, 43] integral equations. Note that in this Fouriertype integral
equation, the scattering matrix stands for a propagator governing the time evolution of the initial spacedependent bright matterwave structure.
To determine conservation laws for the spectral set associated with the inverse scattering problem, we need appropriate compatibility conditions
between equations of the eigensystem Eqs. (7) and (8). For this purpose, we differentiate Eq. (7) with respect to and Eq. (8) with respect
to , and compare the resulting coupled sets. This yields the following secular equation for the spectral parameter
(9) 
It follows that, unlike the homogeneous NLSE for which the IST involves a spectral parameter that is constant in time and space [42, 43], the compatibility equation in Eq. (9) suggests a nonisospectral scattering problem for Eq. (6). An important consequence of the nonuniversality of conservation laws is the possibility of multiple solution to the secular Eq.(9), namely we can peak a solution by the method of separation of variables which amounts to introduce two distinct functions, one depending on space variable and the other on time variable, i.e.
(10) 
In Eq. (9), the latter solution gives rise to the following spacetime evolution equation for the two factors of the spectral parameter
(11) 
It is useful at this step to stress that if was not timedependent,
it would have appeared as a constant factor in Eq. (10) and so
the eigenvalue would be timeindependent but varying with , in
accordance with Eq. (11). Next, one can verify from Eq. (10) that
allowing to change only in time implies that thus the
dependence of the external potential on the space variable should be linear. In this case, Eq. (11) implies that , which is the result of the IST for the linear spacedependent external potential studied in ref. [37].
To find a solution to Eq. (9) in the general case of an external potental wih a linear plus quadratic term, we impose the following constraints
(12) 
where and are two arbitrary real constants. With these constraints Eq. (9) reduces to
(13) 
Hence, the solutions to Eq. (11) are
(14) 
where and are two integration constants. Given that is explicitly defined in (2), equating coefficients of similar terms in this function and in (14) leads to
(15) 
where and are two integration constants. Actually the explicit form for imposes a value on , while we are allowed to set and in this way all coefficients in the spatial component of are unambiguously determined. For the temporal component, we integrate (13) to obtain
(16) 
where
(17) 
with the initial value determined by the scattering potential of the
NIST.
To reconstruct the time evolution of the scattering data, we simply follow the standard considerations of the IST [43], namely, we define a
scattering matrix whose spectral problem involves a finite number () of bound states and a continuum;
(19) 
where the are the amplitudes of the Jost functions [38] which are the solutions to the nonisospectral eigenvalue problems Eqs. (7) and (8), and the indices refer to the discrete eigenstates of the scattering matrix. With the help of the propagator Eq. (LABEL:eq14), we can formulate the inverse transform as follows in terms of the GelfandLevitanMarchenko integral equations [43]
These two integral equations are solved exactly once an explicit shape for the scattering potential is chosen, leading to the following general solution for the nonlinear evolution equation (6)
(21) 
with an arbitrary initial position of the field .
Given that our interest are bright matterwave structures, it is legitimate to choose an initial solution with a pulse shape. Following Satsuma and Yajima [44], such an appropriate initial solution with a pulse shape is given exactly by
(22) 
In quantum mechanics, this sech function represents a reflectionless potential and hence its boundstate spectrum possesses a twofold degenerate mode with the eigenvalues
(23) 
the structure of which suggests that the timedependent part of the nonisospectral eigenvalue of the NIST must be a complex function, namely
(24) 
According to Eq. (23), Eq. (24), and the general expression for given by Eqs. (16)(17), the initial values of the real and imaginary parts of should be
(25) 
Therefore, from (24) we are allowed to write
(26) 
(27) 
Also, the scattering potential (22) permits a complete formulation of the eigenstates of the NIST, including its eigenvalue spectrum and the set of associated eigenfunctions. With this scattering matrix, we construct the timedependent solution to our inhomogeneous NLSE (6) from the integral equation (21). We find
(28) 
with:
(29) 
(30)  
(31)  
where in (29) is a constant amplitude that follows from the explicit form of the scattering matrix via the GelfandLevitanMarchenko integral equations. The exact expression for this parameter resulting from the above treatment is
(32) 
where is the ratio of two arbitrary constant amplitudes of the Jost functions.
3 Variational Versus Exact MatterWave Solitons
In the previous section, the NLSE (1) was solved by following the IST technique, considering a pulseshaped initial profile describing a bright matterwave soliton. Regarding this point, the external potential in the NLSE represents a harmonic hump with an asymmetric shape due to gravity. In practice, such a potential can be associated with a harmonic laser barrier that splits a preformed singlepulse soliton into two pulses with an equal tail, width and hence energy. If one of the two postcreated pulses is allowed to fall freely under gravity, the gravityinduced acceleration will create an extra potential energy, that varies linearly with the pulse position along the direction of free fall. By measuring the energy difference between the nonaccelerated pulse component and the pulse component falling under gravity at some specific position, one can determine exactly the constant of gravity for atoms of a specific species (i.e. of a welldefined mass ) composing a Bose–Einstein condensate. Thus, the problem considered in this work has a direct application in gravity measurements using matterwave interferometry, as explained in more detail in refs. [22, 23].
Although the NIST technique has enabled us find an exact solution to the NLSE (1), in recent studies the problem of bisolitonic wave generation and propagation has most often been considered within the framework of variational approaches [7, 8, 28, 29, 30]. The two most common variational ansatzes used in this context are the supersech and the secondorder Hermite–Gaussian modes.
The supersech mode can be expressed as
(33) 
while the secondorder Hermite–Gaussian mode is given by
(34) 
where , , , and are the initial values (i.e., values at ) of the soliton’s variational parameters, while is the soliton centerofmass position. In the spirit of the variational approach, these two predefined solutions are stationary modes, with (34) representing the second term of a full Hermite–Gaussian expansion of the solution to Eq. (1) without the linear component in the external potential, treating the nonlinear term as a perturbation [27].
In fact, the perturbation theory leading to the Hermite–Gaussian mode (34) is the main factor differentiating the two ansatzes. Indeed, the supersech mode is chosen in such a way as to reproduce the exact NLS singlepulse solution in the absence of external potential. It turns out that the supersech mode has all the features of a true soliton and it is often referred to as a symmetric (or evenparity) mode [7, 8]. In contrast, the Hermite–Gaussian ansatz is the solution to a linear eigenvalue problem with a nonlocal nonlinear spatial external potential. For this reason, the Hermite–Gaussian mode is not a true soliton, despite its localized shape profile. It is now well established that for some values of its variational parameters the Hermite–Gaussian mode is antisymmetric with respect to its argument , which is why it is sometimes also referred to as an oddparity mode [8].
It should be stressed that in our problem, the gravity introduces an extra linear term in the external potential in addition to the quadratic one. Therefore, it is useful to reformulate this variational ansatz. For this purpose we consider a variable change , which transforms the external potential with quadratic plus linear terms given in Eq. (4) into a pure quadratic potential . Doing this, we obtain the same amplitude equation as that solved in ref. [8] for the temporal Hermite–Gaussian mode, which in the present context leads to
(35) 
Concerning the supersech mode, we shall keep the analytical expression (33) since this is actually an arbitrary ansatz as opposed to the Hermite–Gaussian mode which is the solution to an existing equation.
Since our primary goal is to establish a possible consistency between two common variational solutions and the exact soliton solution obtained via the NIST, it will be more practial if we reexpress the exact onesoliton solution given in Eq. (28) as
(36) 
where , a realvalued function representing the NIST bisoliton amplitude, follows straight from Eqs. (29)(32) and thus is given by
(37) 
with
(38) 
the bisoliton centerofmass position
(39) 
its spatial width and the quantity
(40) 
accounting for a residual spatial shift in the bisoliton center. Similarly, the quantity in (37), which follows from Eq. (29) as a prefactor to the sech function, is
(41) 
where
(42)  
(43) 
In the static regime, we set in the above functions, a consideration resulring in two essential remarks:

When we set in Eq. (37), all functions of in this expression become simple parameters. Therefore, we can readily set , , , and , such that the parameters , , , and in the two variational ansatzes are equal to the equivalent parameters in the exact onesoliton solution in the static regime.

According to Eq. (26), the quantity is always positive for finite values of . It follows that the parameter , obtained from Eq. (40) with , can take positive and negative values depending on the value of given in formula (32). Quite remarkably, while this sign change preserves the bisoliton profile of the exact soliton solution, numerical analysis of the supersech mode reveals that the change in from positive to negative values gives rise to either a single pulse or a bisoliton. Therefore, we can conclude that in general the bisoliton profile of the propagating supersech mode is determined by appropriate values of its variational parameter , obtained by solving the associated variational equation.
We start our analysis by comparing intensity profiles of the supersech mode with those of the exact bisoliton. The three left graphs of figure 1 represent the amplitudes of the exact bisoliton solution, while the three right graphs are the amplitudes of the supersech ansatz for and and for three different values of , i.e., , , and . As one can see, while the amplitude of the NIST solution always has a bisolitonic shape, for positive the supersech mode has a singlepulse shape with an increasingly broad peak as increases.
Oppositely, when is decreased in the negative branch, the singlepulse shape
breaks up into a bisoliton, as is apparent in the bottom left
graph of figure 1. The emerging twin pulses are increasingly separated with decreasing
in the negative branch. Note that in both cases, the peak intensities of the twin pulses are equal.
In addition to their equal peak intensities, the two pulses seen in figure 1 are always symmetric with respect to the fixed center of mass . This
feature, together with the equal peak intensities, confers these specific bisolitonic structures with a unique character, which we term the ”genuine pulse twinning” state. The main distinctive feature is that the bisoliton energy is always exactly twice the energy of each of the two constituent
pulses, even in the regime where the twin pulses partially overlap.
In our next analysis, we compare the intensity profile of the exact NIST bisoliton with
the intensity profiles of the Hermite–Gaussian and supersech modes. The three left graphs of figures 2 and 3 are intensity profiles
of the NIST soliton (solid curves) and the
supersech mode (dashed curves), while the three right graphs are intensity profiles of
the NIST soliton (solid curves) and the Hermite–Gaussian mode (dashed
curves), both plotted in the same graph. The curves in figure 2 are for and , while the
curves in figure 3 are for and a relatively larger value of (i.e., ). Values of
are selected in a range where both the supersech and Hermite–Gaussian modes display
bisoliton shapes.
According to the left graphs of figures. 2 and 3, when is decreased in the
negative branch, the NIST soliton and the supersech mode tend to two qualitatively identical bisolitonic structures. In contrast to this behavior, the right graphs suggest a tendency to
identical bisoliton profiles provided is increased in the positive
branch. Thus, the ranges of values in which the
evenparity and oddparity bisolitons are both consistent with the shape profiles of the exact NIST soliton do not match.
In the above analysis, we assumed corresponding to the context of beam splitting without postgravitational acceleration of the output matterwave twin pulses. Now we turn to the case that atoms experience gravity. For this case we have shown above that it is possible to mathematically account for the effect of gravity on the analytical expression for the Hermite–Gaussian mode (see Eq. 35), but this is not possible for the supersech ansatz, which is actually an arbitrarily chosen variational ansatz. Therefore, in this second analysis, we focus on the comparision of the Hermite–Gaussian mode, now given by Eq. (35), with the exact onesoliton solution. In figure 4 we show the intensity profile of the Hermite–Gaussian mode for with , (two left graphs), and , (two right graphs) in order to see the effects of the linear term in the external potential on the Hermite–Gaussian quasibisoliton shape.
For the exact bisoliton, in figures 5 and 6 we plotted its intensity profiles for , , and , , respectively, for three distinct negative (left graphs) and positive (right graphs) values of .
When both and are relatively small (Fig. 5), the bisoliton feature survives but the two pulses have different widths and tails. In addition, the center positions of the twinned pulses are no longer symmetric with respect to , as was the case in the absence of a gravity potential. Instead we have two asymmetric pulses. On the other hand, when the linear coefficient is relatively large compared with the quadratic coefficient (Fig. 6), the bisoliton shape is destroyed for certain values of . Given that for the exact bisoliton, depends mainly on the characteristic parameters of the matterwave evolution equation, we can exqctly determine the existence condition for a matterwave bisoliton structure. These features of the exact bisoliton solution, namely, the shift of pulse center positions with varying characteristic parameters of the system, cannot be accounted for by the Hermite–Gaussian mode as is apparent in figure 4. Indeed, figure 4 clearly shows that although the intensities of the two pulses are not equal, their amplitudes persistently overlap for all values of .
4 Numerical Simulations
In the previous section, an extensive analysis of the consistency between the exact soliton solution to the GP Eq. (6), on one hand, and two variational ansatzes to the same equation, on the other hand, was carried out. We found that the two variational fields agree qualitatively with the exact solution only in restricted ranges of their characteristic parameters. Most importantly, the parameter values for which the Hermite–Gaussian and the supersech modes display the required doublepulse profile are distinct. However, the exact solution is always a doublepulse structure, with the pulse center positions and intensities changing with the linear and quadratic coefficients of the external potential.
Nevertheless, the analysis in the previous section considered only static profiles of the fields. This is because the dynamical properties of the two variational fields can be determined only by solving appropriate variational equations, which is beyond the scope of the present work. In fact, our main interest was the consistency of the shape profiles of the three solutions. Concerning the system dynamics, the most relevant information about the behaviors of the matterwave field is expected from the exact solution to the equation of motion, i.e., the GP equation. In this respect, we simulated the GP equation (6) by following a splitstep scheme in the time domain, assuming that the spatial profile of the field is given by Eq. (36) at . With this consideration, we can vary the characteristic parameter , in addition to the linear and quadratic coefficients of the external position.
Figures 7 (), 8 (), and 9 () display profiles of the field intensity in spacetime coordinates for some sets of values of characteristic model parameters.
Figures 7, which corresponds to the numerical results for , shows that in the absence of gravity, the matterwave field has a symmetric doublepulse structure throughout its propagation time. When (i.e., in the absence of gravity), the two copropagating pulses form a twin pair with the same intensity, but the separation between their peaks seems to change only weakly. The same behavior emerges in figure 8 for a nonzero but positive value of . However, when is negative, the peak separation manifestly varies in time. More explicitely, the two pulses move away from each other with time, while the heights of their peaks increase. Note that the quantity for the exact soliton solution to the GP equation is defined in Eq. (40) it is indeed negative by virtue of the first term in the righthand side of Eq. (40) if is very small compared with .
5 Discussion and Concluding Remarks
Matterwave interferometry is one of the most active current fields of investigation in quantum metrology, because of its great advantage related to the finite mass of the field constituents that permits the probing of quantum properties of matter, and force fields such as gravity, mattermatter interactions and so forth.
In a doubleslit experiment involving a matterwave field to measure the gravity force, for instance, one can assume a single BEC initially prepared in an optical barrier such as a faroff resonant laser barrier [22, 23], which splits the single BEC into two separate qualitatively and quantitatively identical pulsed fields. When one of the two pulses is released under gravity, this causes asymmetry of the effective optical potential consisting of the antiharmonic optical trap and the gravitational field, and the twopulse matterwave soliton forms two copropagating unbalanced intensity pulses as a result of gravitational acceleration of one of the two pulses. By measuring the potential energy difference between the two copropagating pulses, one can, for instance, estimate the kinetic energy due to this gravitational acceleration. In this regard, the analysis of sect. 3 shows that when the coefficient of the linear term in the potential is zero, the two pulses are perfectly twinned, indicating the absence of gravitational acceleration of pulses. When is nonzero, the two pulses become increasingly unbalanced as is increased. Noting that , according to Eq. (5), is proportional to the constant of gravity and the mass of atom species in the condensate, the behavior of the accelerated pulse fraction under the variation of reflects an enhancement of the gravitational effect as atoms forming the condensate are heavier.
In modeling the doubleslit interferometry, an appropriate description of the splitting of a single BEC into two copropagating matterwave pulses, taking into account the necessity to accurately control the spatial and temporal locations of each pulse, is essential in probing spatially and temporally varying forces or interactions. While a doublepulse matterwave field has usually been approximated in variational treatments either by the Hermite–Gaussian wave or supersech mode, neither of these two fields is an exact solution to the GP equation governing the system dynamics.
In this work, we addressed the issue of the consistency of these variational fields versus an exact onesoliton solution to the GP equation obtained by means of an inversescattering transform. The following general picture emerges from our analysis:

The supersech mode has a bisoliton shape with respect to the spatial coordinate when a phase factor in its argument is negative. This bisoliton profile becomes sharper, and pulses in the bound state increasingly well seperated, as is decreased in the negative branch. In other words, when is negative but close to zero, the supersech mode consists of two partially overlapping pulses the separation of which increases as decreases in the negative branch.

The Hermite–Gaussian mode has a bisoliton profile for all nonzero values of its characteristic parameters. However, unlike the supersech mode, the Hermite–Gaussian bisoliton is odd with respect to the spatial coordinate. The Hermite–Gaussian bisoliton describes two pulses that are always partially embedded in each other in the static regime. Therefore, only by propagation can the two pulses deploy themselves into two more or less independent individual entities, as has been observed in numerical simulations in the context of dispersionmanaged optical fibers [8, 45].

The NIST soliton solution possesses a permanent bisolitonic profile, in addition, its propagation leads to two individual pulseshaped entities moving away from each other with different intensities.
The results of numerical simulations have emphasized the virtues of the exact soliton solution, and particularly its ability to provide a coherent description of the formation and propagation of two coexisting pulses with distinct dynamical parameters. We therefore conclude that although variational ansatzes such as the Hermite–Gaussian and supersech modes are relevant to a qualitative description of atom interferometry involving twofield matterwave modes, an exact solution to the governing GP equation including a beam splitter and gravity should provide the most physically relevant insight into the system dynamics.
The authors thank the referee for drawing their attention to the important work of Satsuma and Yajima (ref. 44) on initialvalue problems for the nonlinear Schrödinger equation. This work was partially supported by the Academic of Science for the Developing World (TWAS). The authors are grateful to the Physical Society of Japan for financial support in publication.
References
 D. Y. Tang, W. S. Man, H. Y. Tam and P. D. Drummond: Phys. Rev. A 64 (2001) 033814.
 D. Y. Tang, L. M. Zhao, B. Zhao and A. Q. Liu: Phys. Rev. A 72 (2005) 043816.
 D. Y. Tang, L. M. Zhao and B. Zhao: Appl. Phys. B 80 (2005) 239.
 D. Y. Tang, B. Zhao, D. Y. Shen, C. Lu, W. S. Man and H. Y. Tam: Phys. Rev. A 66 (2002) 033806.
 J. M. SotoCrespo, N. Akhmediev, Ph. Grelu and F. Belhache: Opt. Lett. 28 (2003) 1757.
 D. Jr. Fandio Jubgang, A. M. Dikandé and A. SundaMeya: Phys. Rev. A 92 (2015) 053850.
 S. Kumar and A. Hasegawa: Opt. Lett. 22 (1997) 372.
 C. Paré and P. A Boulanger: Opt. Commun. 168 (1999) 103.
 A. Maruta, T. Inoue, Y. Nonaka and Y. Yoshika: Electron. Lett. 37 (2001) 1357.
 A. Maruta, T. Inoue, Y. Nonaka and Y. Yoshika: IEEE J. Sel. Top. Quantum Electron. 8 (2002) 640.
 B. A. Malomed: Opt. Commun. 136 (1997) 313.
 R. K. Jackson, Christopher K.R.T. Jones and V. Zharnitsy: Physica D 190 (2004) 63.
 BaoFeng Fang and Boris A. Malomed: Opt. Commun. 229 (2004) 173.
 P. Kockaert and S. Larochelle: Opt. Commun. 246 (2005) 185.
 M. Stratmann, T. Pagel and F. Mitschke: Phys. Rev. Lett. 95 (2005) 143902.
 A. Chong, W. H. Renninget and F. W. Wise: Opt. Lett. 33 (2008) 1717.
 E. Infeld, M. Matuszewski, C. Shino and M. Trippenbach: Optica Applicata 36 (2006) 575.
 E. Infeld, M. Matuszewski and M. Trippenbach: J. Phys. B 39 (2006) L113.
 H. Fujishima, M. Okumura, M. Mine and T. Yajima: J. Phys. Soc. Jpn. 81 (2012) 104003.
 Y. Shin, M. Saba, T.A. Pasquini, W. Ketterle, D.E. Pritchard, and A. E. Leanhardt: Phys. Rev. Lett. 92 (2004) 050405.
 Y. Shin, C. Sanner, G. B. Jo, T.A. Pasquini, M. Saba, W. Ketterle, D. E. Pritchard, M. Vengalattore and M. Pentiss: Phys. Rev. A 72 (2005) 021604.
 L. Pezzé, L. A. Collins, A. Smerzi, G. P. Berman and A. R. Bishop: Phys. Rev. A 72 (2005) 043612.
 L. Pezzé, A. Smerzi, G. P. Berman, A. R. Bishop and L. A. Collins: Phys. Rev. A 74 (2006) 033610.
 D. M. StamperKurn, M. R. Andrews, A. P. Chikkatur, S. Inouye, H.J. Miesner, J. Stenger and W. Ketterle: Phys. Rev. Lett. 80 (1998) 2027.
 S. Raghavan, A. Smerzi, S. Fantoni and S. R. Shenoy: Phys. Rev. A 59 (1999) 620.
 S. K. Turitsyn, B. G. Bale and M. P. Fedoruk: Phys. Rep. 521 (2012) 135.
 T. I. Lakoba and D. J. Kaup: Phys. Rev. E 58 (1998) 6728.
 A. Klein, D. Jaksch, Y. Zhang and W. Bao: Phys. Rev. A 76 (2007) 043602.
 D. A. Zezyulin, G. L. Alfimov, V. V. Konotop and V. M. PérezGarcía: Phys. Rev. A 78 (2008) 013606.
 T. P. Billam, S. A. Wrathmall and S. A. Gardiner: Phys. Rev. A 85 (2012) 013627.
 Y. S. Kivshar and B. A. Malomed: Rev. Mod. Phys. 61 (1989) 763.
 A. M. Dikandé: J. Math. Phys. 49 (2008) 073520.
 M. Wadati: J. Phys. Soc. Jpn. 68 (1999) 2543.
 M. Wadati: J. Phys. Soc. Jpn. 70 (2001) 60.
 Y. Le Coq, J. H. Thywissen, S. A. Rangwala, F. Gerbier, S. Richard, G. Delannoy, P. Bouyer and A. Aspect: Phys. Rev. Lett. 87 (2001) 170403.
 N. P. Robins, A. K. Morrison, J. J. Hope and J. D. Close: Phys. Rev. A 72 (2005) 631606.
 H. H. Chen and C. S. Liu: Phys. Rev. Lett. 37 (1976) 693.
 R. Balakrishnan: Phys. Rev. A 32 (1985) 1144.
 Z. Xu, L. Li, Z. Li, G. Zhou and K. Nakkeeran: Phys. Rev. E 68 (2003) 046605.
 C. C. Mak, W. Chow and K. Nakkeeran: J. Phys. Soc. Jpn. 74 (2005) 1449.
 K. Nakkeeran: J. Phys. A 34 (2001) 5111.
 M. J. Ablowitz and H. Segur: Solitons and the Inverse Scattering Transform (SIAM, Philadelphia, PA, 1981).
 R. K. Dodd, J. C. Eilbeck, J. D. Gibbon and H. C. Morris: Solitons and Nonlinear Wave Equations (Academic, New York, 1982).
 J. Satsuma and N. Yajima: Suppl. Prog. Theor. Phys. 55 (1974) 284.
 C. Paré: Opt. Commun. 154 (1998) 9.