Matter-Wave Fields for Double-Slit Atom Interferometry: Variational Versus Exact Solitons\abst
A major challenge in the theoretical modeling of double-slit interferometry involving matter-wave 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 single-hump 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 super-sech modes. From numerical simulations, evidence is given of the exact solution as being the most appropriate matter-wave structure that provides a coherent description of the generation and spatio-temporal evolution of matter-wave optical fields in a hypothetical implementation of double-slit atom interferometry.
Bound-soliton pairs provide interesting transmission channels for multiplexed high-intensity 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 single-pulse 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 matter-wave pulse, which creates two co-propagating matter-wave pulses.
Bisolitons are of fundamental importance in BECs since they provide appropriate wave structures needed for the experimental achievement of an atomic equivalent of double-slit holography [22, 23]. In such an experiment, we can imagine two spatially and temporally entangled matter-wave pulses generated from the splitting of a single condensate into two atomic populations, namely in nearly degenerate hyperfine levels (see for example ref. ). Generally a far-off resonant laser barrier provides an ideal scheme for such splitting [22, 23]. In theory, this most often has been described by considering a single matter-wave 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. , in terms of a single macroscopic atomic population moving in a double-well 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. , 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 matter-wave 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 super-sech 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 .
In this work, we address the issue by considering two distinct variational solutions, namely, the super-sech and Hermite–Gaussian modes [19, 28, 29, 30] on one hand, and the exact one-soliton solution to the associated GP equation obtained by a non-isospectral inverse-scattering transform (NIST) method  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 post-created bisoliton with respect to the other due to the free-fall 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 matter-wave interferometry as previously emphasized (see for example refs. [22, 23]).
2 Exact One-Soliton 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 mean-field 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
where is the mean–field interatomic interaction coefficient related to the s-wave scattering length. is the total external potential for the nonlinear Schrödinger-type Eq. (1), it is given as the sum of the repulsive harmonic potential and the linear potential accounting for the gravitational field
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
With these new variables, the GP equation in Eq.(1) becomes
which is simply the self-focusing NLSE with a nonlocal (i.e., spatially varying) external potential . When and , Eq. (6) describes the macroscopic dynamics of a free-falling BEC. The same equation was proposed to model the dynamics of plasmas in a gravitational field, and its exact one and n-soliton solutions have been obtained  by means of the Inverse-Scattering Transform (IST). On the other hand, when but , Eq. (6) becomes the NLSE with an external space-varying repulsive harmonic potential. Exact soliton solutions for this latter physical context have also been obtained, using an improved IST technique that overcomes the non-isospectral character of the associated Lax-pair-type eigenvalue problem .
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 . 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 
where and are space-time-dependent functions that form a two-component 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 Lax-pair 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., pulse-shaped) 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  of the initial scattering potential , and from this choice the time evolution of the scattering data can easily be computed via Gelfand-Levitan-Marchenko [42, 43] integral equations. Note that in this Fourier-type integral equation, the scattering matrix stands for a propagator governing the time evolution of the initial space-dependent bright matter-wave 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
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 non-isospectral scattering problem for Eq. (6). An important consequence of the non-universality 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.
In Eq. (9), the latter solution gives rise to the following space-time evolution equation for the two factors of the spectral parameter
It is useful at this step to stress that if was not time-dependent,
it would have appeared as a constant factor in Eq. (10) and so
the eigenvalue would be time-independent 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 space-dependent external potential studied in ref. .
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
where and are two arbitrary real constants. With these constraints Eq. (9) reduces to
Hence, the solutions to Eq. (11) are
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
with the initial value determined by the scattering potential of the
To reconstruct the time evolution of the scattering data, we simply follow the standard considerations of the IST , namely, we define a scattering matrix whose spectral problem involves a finite number () of bound states and a continuum;
where the are the amplitudes of the Jost functions  which are the solutions to the non-isospectral 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 Gelfand-Levitan-Marchenko integral equations 
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)
with an arbitrary initial position of the field .
Given that our interest are bright matter-wave structures, it is legitimate to choose an initial solution with a pulse shape. Following Satsuma and Yajima , such an appropriate initial solution with a pulse shape is given exactly by
In quantum mechanics, this sech function represents a reflectionless potential and hence its boundstate spectrum possesses a twofold degenerate mode with the eigenvalues
the structure of which suggests that the time-dependent part of the non-isospectral eigenvalue of the NIST must be a complex function, namely
Therefore, from (24) we are allowed to write
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 time-dependent solution to our inhomogeneous NLSE (6) from the integral equation (21). We find
where in (29) is a constant amplitude that follows from the explicit form of the scattering matrix via the Gelfand-Levitan-Marchenko integral equations. The exact expression for this parameter resulting from the above treatment is
where is the ratio of two arbitrary constant amplitudes of the Jost functions.
3 Variational Versus Exact Matter-Wave Solitons
In the previous section, the NLSE (1) was solved by following the IST technique, considering a pulse-shaped initial profile describing a bright matter-wave 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 single-pulse soliton into two pulses with an equal tail, width and hence energy. If one of the two post-created pulses is allowed to fall freely under gravity, the gravity-induced 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 non-accelerated 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 well-defined mass ) composing a Bose–Einstein condensate. Thus, the problem considered in this work has a direct application in gravity measurements using matter-wave 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 super-sech and the second-order Hermite–Gaussian modes.
The super-sech mode can be expressed as
while the second-order Hermite–Gaussian mode is given by
where , , , and are the initial values (i.e., values at ) of the soliton’s variational parameters, while is the soliton center-of-mass 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 .
In fact, the perturbation theory leading to the Hermite–Gaussian mode (34) is the main factor differentiating the two ansatzes. Indeed, the super-sech mode is chosen in such a way as to reproduce the exact NLS single-pulse solution in the absence of external potential. It turns out that the super-sech mode has all the features of a true soliton and it is often referred to as a symmetric (or even-parity) 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 odd-parity mode .
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.  for the temporal Hermite–Gaussian mode, which in the present context leads to
Concerning the super-sech 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 one-soliton solution given in Eq. (28) as
the bisoliton center-of-mass position
its spatial width and the quantity
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 one-soliton 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 super-sech 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 super-sech 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 super-sech 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 super-sech 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 super-sech mode has a single-pulse shape with an increasingly broad peak as increases.
Oppositely, when is decreased in the negative branch, the single-pulse 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 super-sech modes. The three left graphs of figures 2 and 3 are intensity profiles of the NIST soliton (solid curves) and the super-sech 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 super-sech 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 super-sech 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
even-parity and odd-parity 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 post-gravitational acceleration of the output matter-wave 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 super-sech 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 one-soliton 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 quasi-bisoliton shape.
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 matter-wave evolution equation, we can exqctly determine the existence condition for a matter-wave 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 super-sech modes display the required double-pulse profile are distinct. However, the exact solution is always a double-pulse 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 matter-wave 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 split-step 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 space-time 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 matter-wave field has a symmetric double-pulse structure throughout its propagation time. When (i.e., in the absence of gravity), the two co-propagating 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 right-hand side of Eq. (40) if is very small compared with .
5 Discussion and Concluding Remarks
Matter-wave 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, matter-matter interactions and so forth.
In a double-slit experiment involving a matter-wave field to measure the gravity force, for instance, one can assume a single BEC initially prepared in an optical barrier such as a far-off 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 two-pulse matter-wave soliton forms two co-propagating unbalanced intensity pulses as a result of gravitational acceleration of one of the two pulses. By measuring the potential energy difference between the two co-propagating 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 double-slit interferometry, an appropriate description of the splitting of a single BEC into two co-propagating matter-wave 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 double-pulse matter-wave field has usually been approximated in variational treatments either by the Hermite–Gaussian wave or super-sech 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 one-soliton solution to the GP equation obtained by means of an inverse-scattering transform. The following general picture emerges from our analysis:
The super-sech 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 super-sech 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 super-sech 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 dispersion-managed optical fibers [8, 45].
The NIST soliton solution possesses a permanent bisolitonic profile, in addition, its propagation leads to two individual pulse-shaped 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 co-existing pulses with distinct dynamical parameters. We therefore conclude that although variational ansatzes such as the Hermite–Gaussian and super-sech modes are relevant to a qualitative description of atom interferometry involving two-field matter-wave 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 initial-value 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.
-  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. Soto-Crespo, N. Akhmediev, Ph. Grelu and F. Belhache: Opt. Lett. 28 (2003) 1757.
-  D. Jr. Fandio Jubgang, A. M. Dikandé and A. Sunda-Meya: 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.
-  Bao-Feng 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. Stamper-Kurn, 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érez-Garcí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.