Fully three dimensional breather solitons can be created using Feshbach resonance
M.~Matuszewski
We investigate the stability properties of breather solitons in a three-dimensional Bose-Einstein Condensate with Feshbach Resonance Management of the scattering length and confined only by a one dimensional optical lattice. We compare regions of stability in parameter space obtained from a fully 3D analysis with those from a quasi two-dimensional treatment. For moderate confinement we discover a new island of stability in the 3D case, not present in the quasi 2D treatment. Stable solutions from this region have nontrivial dynamics in the lattice direction, hence they describe fully 3D breather solitons. We demonstrate these solutions in direct numerical simulations and outline a possible way of creating robust 3D solitons in experiments in a Bose Einstein Condensate in a one-dimensional lattice. We point other possible applications.
FRM_let.tex
\documentclass[aps,prl,twocolumn,superscriptaddress,showpacs]{revtex4}
\usepackage{graphicx}
\bibliographystyle{apsrev}

\newcommand{\rbm}[1]{\bm{\mathrm{#1}}}

\begin{document}

\title{Fully three dimensional breather solitons can be created using Feshbach resonance}

\author{M.~Matuszewski}
\affiliation{Institute of Theoretical Physics, Physics Department, Warsaw University,
Ho\.{z}a 69, PL-00-681 Warsaw, Poland}

\author{E.~Infeld}
\affiliation{Soltan Institute for Nuclear Studies, Ho\.{z}a 69,
PL-00-681 Warsaw, Poland}

\author{B.~A.~Malomed}
\affiliation{Department of Interdisciplinary Sciences, School of Electrical
Engineering, Faculty of Engineering, Tel Aviv University, Tel Aviv 69978, Israel}

\author{M.~Trippenbach}
\affiliation{Institute of Theoretical Physics, Physics Department, Warsaw University,
Ho\.{z}a 69, PL-00-681 Warsaw, Poland}

\begin{abstract}
We investigate the stability properties of breather solitons in a three-dimensional
Bose-Einstein Condensate with Feshbach Resonance Management of the scattering length and
confined only by a one dimensional optical lattice. We compare regions of stability in
parameter space obtained from a fully 3D analysis with those from a quasi two-dimensional
treatment. For moderate confinement we discover a new island of stability in the 3D case,
not present in the quasi 2D treatment. Stable solutions from this region have nontrivial
dynamics in the lattice direction, hence they describe fully 3D breather solitons. We
demonstrate these solutions in direct numerical simulations and outline a possible way of
creating robust 3D solitons in experiments in a Bose Einstein Condensate in a
one-dimensional lattice. We point other possible applications.
\end{abstract}
\pacs{03.75.Lm, 05.45.Yv}

\maketitle

The creation of Bose--Einstein condensates (BEC) in vapors of alkali metals has opened an
excellent opportunity to investigate nonlinear interactions of atomic matter waves. An
important challenge of practical interest is to develop methods to create and control
matter-wave solitons. One dimensional dark \cite{1fsol}, bright \cite{2bright}, and
gap-mode \cite{ober} solitons have already been observed. A promising approach to obtain
multidimensional solitons consists in varying the scattering length {\it a} of
interatomic collisions. This can be achieved by means of sweeping an external magnetic
field through the $a=0$ point (the point, where the scattering length vanishes). This
point occurs close to the Feshbach resonance \cite{3Inouye}. The application of an ac
magnetic field may induce a periodic modulation of {\it a}, opening a way to
``Feshbach-resonance management" (FRM) \cite{4Greece}. A noteworthy FRM-induced effect is
the possibility of creating self-trapped oscillating BEC solitons (breathers) without an
external trap in the 2D case.  The underlying mechanism is fast modulations creating an
effective potential on a slower timescale. This potential can stabilize the soliton. The
BEC model based on the Gross-Pitaevskii equation (GPE) with harmonic modulation of {\it
a} was investigated in Refs.~\cite{6Saito,8Abdul,9VPG}. The conclusion was that FRM
renders it possible to stabilize 2D breather solitons even without the use of an external
trap. According to these references, 3D breathers require at least a tight, one
dimensional harmonic trap \cite{Gerlitz}, practically reducing the problem to 2D. Later
on we will call this approach a quasi two dimensional (Q2D) treatment. It will be defined
more precisely below.

Recently \cite{10EL} we demonstrated  that quasi 2D solitons can be stabilized by a
combination of FRM and a strong 1D optical lattice (1D OL), instead of a 1D harmonic trap
\cite{6Saito,8Abdul,9VPG}. By a ``strong lattice" we mean one in which the atoms in
neighboring cells cannot interact. This issue has practical relevance, as a 1D OL can
easily be created, illuminating the BEC by a pair of counterpropagating laser beams that
form a periodic interference pattern \cite{11OL}. The lattice will be
weak or strong depending on the intensity of the laser light. In fact, it is easier to realize a
tight confinement configuration in an optical lattice than in an harmonic trap. Hence
this environment may be more friendly for creating quasi 2D solitons.

In this Letter, we demonstrate that the combined OL-FRM stabilization of 3D solitons is possible
even in a weak lattice, when atoms confined in different cells interact. By analyzing the
stability charts in configuration space we discover two distinct regions where stable solutions
exist. The first of these regions has its counterpart in the Q2D treatment. The other region
appears when the frequency of modulation exceeds a critical value dependent on the
confining potential.
%the lowest excitation frequency of the confining potential.
It is not present in the Q2D treatment; hence it corresponds to fully 3D solitons. In the
limit of tight confinement the latter region moves to extremely high frequencies and the
Q2D stability chart is recovered.

We describe our system by the GPE in physical units, including a
time-dependent (FRM-controlled) scattering length $a(\tau)$ and an
external potential $\tilde{U}(\mathbf{r},\tau)$
\begin{equation}
i \hbar\frac{\partial \Psi}{\partial \tau}=\left[ -\frac{\hbar^2}{2m}\nabla
^{2}+\tilde{U}(\mathbf{r},\tau)+ \frac{4 \pi a(\tau) \hbar^2}{m}|\Psi|^{2}\right] \Psi .
\label{NLS0}
\end{equation}
Initially the BEC is in the ground state of a radial (2D) parabolic trap with frequency
$\tilde{\omega}_{\perp}$, supplemented, in the longitudinal direction, by ``end caps"
induced by transverse light sheets. The configuration is much like the one used to create
soliton trains in a Li${^{7}}$ condensate \cite{2bright}. A 1D lattice potential in the
axial direction is adiabatically turned on from $\tilde{\varepsilon}=0$ to
$\tilde{\varepsilon}=\tilde{\varepsilon}_{\mathrm{f}}$,  see Fig.~\ref{fig1}. Thus, the full
potential is
\begin{eqnarray}
\tilde{U}(\mathbf{r},\tau)&=&\tilde{\varepsilon} (\tau)\left[ 1-\cos (2 \pi z/\lambda)\right] + \nonumber \\
&+&f(\tau)\left[(m/2)\tilde{\omega}_{\perp}^{2}\varrho ^{2}+\tilde{U}_{0}(z)\right],
\label{V0}
\end{eqnarray}
where $\lambda$ is the lattice spacing, $\varrho$ is the radial variable in the plane
transverse to $z$, and the axial ``end-cap'' potential, $\tilde{U}_{0}(z)$ is
approximated by a sufficiently deep one dimensional rectangular potential well. The width
of the well determines the number of peaks in the finally established structure. The
$f(\tau)$ is a switching-off function (see Fig.~\ref{fig1}). We introduce dimensionless
variables $\mathbf{x}=(\pi/\lambda)\mathbf{r}$, $t=\tau \omega$, $\psi=\Psi
\sqrt{\lambda^3/(N\pi^3)}$, where $\omega=\pi^2 \hbar/(m \lambda^2)$ and $N$ is the
number of atoms. We also define $g=4 \pi^2 N a /\lambda$,
$\varepsilon=\tilde{\varepsilon}/(\omega\hbar)$, $\omega_{\perp}=\tilde{\omega}_{\perp
}/\omega$, and $U_0=\tilde{U}_{0}/(\omega\hbar)$ to obtain
\begin{equation}
i\frac{\partial \psi}{\partial t}=\left[ -\frac{\nabla ^{2}}{2}+U(\mathbf{x},t)+g(t)|\psi
|^{2}\right] \psi , \label{NLS}
\end{equation}
and the potential in the form
\begin{equation}
U(\mathbf{x},t)=\varepsilon (t)\left[ 1-\cos (2z)\right] +f(t)\left[\frac{1}{2}\omega _{\perp
}^{2}\varrho ^{2}+U_{0}(z)\right].  \label{V}
\end{equation}
The nonlinear interaction coupling is described by
\begin{equation}
g(t)=g_{0}(t)+g_{1}(t)\sin (\Omega t),  \label{g}
\end{equation}
and the dimensionless modulation frequency by
$\Omega=\tilde{\Omega}/\omega$. Initially $g_1(0)=0$ and $g(0)=g_0(0)>0$.
At some moment $t_{1}$, we begin to linearly decrease $g_0(t)$. It
vanishes at time $t_{2}$, and remains zero up to $t_{3}$, when we
start to gradually switch on the rapid FRM modulation of {\it a}. In
the interval $[t_3,t_4]$, $g_0(t)$ decreases linearly from zero to a
negative $g_{0 \mathrm{f}}$ and the amplitude of the modulation
$g_1(t)$ increases from zero to $g_{1 \mathrm{f}}$, see
Fig.~\ref{fig1}. Simultaneously,  both the radial confinement and
end-caps are gradually switched off.
%(see the the behavior of the $f(t)$ function in Fig. \ref{fig1})
At times $t>t_{4}$, $g(t)$ oscillates with a constant amplitude $g_{1 \mathrm{f}}$ around a
negative average value $g_{0 \mathrm{f}}$. Consequently, a soliton so created, if any,
is supported by the combination of the 1D lattice and FRM.

\begin{figure}[tbp]
\includegraphics[width=8.5cm]{schemeColor.eps}
\caption{(Color online). The time dependence of the nonlinear coefficient, $g$, switching-off function, $f(t)$,
and optical-lattice strength, $\protect\varepsilon ,$ in the numerical experiment which leads to
the establishment of stable 3D breathing solitons supported by the combination of the quasi-1D
lattice and Feshbach-resonance management (FRM). The shaded area indicates rapid oscillations of
$g(t)$, which account for the FRM. \label{fig1}}
\end{figure}

Numerical experiments following the path outlined in Fig.~\ref{fig1} indicate that it is possible
to create stable solitons \cite{10EL} (see inset to Fig.\ref{fig2}). Before showing the
results, we first resort to the variational approximation (VA) in order to predict conditions on
the modulation frequency and the size of the negative average nonlinear coefficient $g_{0 \mathrm{f}}$,
necessary to support 3D solitons.

%Note that our system has an alternative physical interpretation in terms of nonlinear optics
%\cite{5Isaac,10EL}.

The VA can be applied to the description of BEC dynamics under diverse circumstances
\cite{6Saito,8Abdul,9VPG,10EL,12BBB,Progress}. Equation (\ref{NLS}) is derived from the Lagrangian
density
\begin{eqnarray}
{\cal{L}}= i(\psi _{t}^{\ast }\psi -\psi _{t}^{\ast }\psi )-|\psi _{\varrho }|^{2}-|\psi _{z}|^{2}
-g(t)|\psi |^{4} - 2U|\psi |^{2} \label{L}
\end{eqnarray}
We use VA for $t>t_3$ and choose a complex Gaussian ansatz for the solution for one
lattice cell. The amplitude is $A(t)$, radial and axial widths are $W(t)\ $ and $V(t)$
respectively, and $b(t)$ and $\beta (t)$ are the corresponding chirps
\begin{equation}
\psi (\mathbf{r},t)=Ae^{\left[ -\varrho^{2}(1/2W^{2}+ib) -z^{2}( 1/2V^{2}+i\beta)+i\phi\right]}.
\label{ansatz}
\end{equation}
The reduced Lagrangian can be found upon substituting (\ref{ansatz}) into (\ref{L}) and
integrating over space.
%% \begin{eqnarray}
%% L & = & A^2W^2V \, \bigg[ \frac{ A^2 g}{4\sqrt{2}}+\frac{\varepsilon}{ e^{V^2}} -
%% \frac{f(t)\omega_\perp^2W^2}{2} -\frac{1}{4 V^2} +  \nonumber  \\ & - & \frac{V^2 }{4}
%% \left(\beta^2 + \dot{\beta}\right) -\frac{1}{2 W^2}-\frac{W^2 }{2}\left(b^2+\dot{b} \right)-
%% \dot{\phi} \bigg]
%% \end{eqnarray}
By varying this reduced Lagrangian with respect to $\phi$ we obtain the constant
$E=A^{2}W^{2}V =\pi ^{-3/2}\int_{cell}|\psi |^{2}d\mathbf{r}=\pi ^{-3/2}/n$, where the
integral extends over one cell of the lattice and $n$ is the number of occupied lattice
cells. Notice that the total number of atoms is included in the definition of the
nonlinear coupling $g(t)$ and the total wavefunction is normalized to unity. When the
other four variational equations are derived, we can deduce two dynamical equations for
the widths:
\begin{eqnarray}
\ddot{W} &=&\frac{1}{W^{3}}-f(t)\omega _{\perp }^{2} W+ \frac{Eg(t)}{\sqrt{8}W^{3}V},
\label{variat1} \\
\ddot{V} &=&\frac{1}{V^{3}}-4\varepsilon_{\mathrm{f}} V\exp \left( -V^{2}\right)
+\frac{Eg(t)}{\sqrt{8}W^{2}V^{2}}.  \label{variat2}
\end{eqnarray}
These equations describe the dynamics of a single peak, and with one term slightly altered can
be applied to the problem of a BEC confined in a 1D harmonic trap \cite{6Saito,8Abdul,9VPG}.

In the corresponding Q2D treatment we drop the $z$ dimension in Eq.~(\ref{NLS}). It is
assumed that in this direction the profile of the wavefunction is fixed and reproduces
the ground state $\psi_0$ of the single lattice cell (Wannier function) or harmonic
potential, as in Refs.~\cite{6Saito,9VPG}. The reduced potential in 2D GP will take the
form $U(\varrho,t)= f(t)(1/2)\omega _{\perp }^{2}\varrho ^{2}$. In VA we take $V \equiv
V_0$ from $4\varepsilon_{\mathrm{f}} V_{0}^{4}\exp \left( -V_{0}^{2}\right)=1$, and only
solve equation (\ref{variat1}). In numerical simulations we rescale the nonlinear
coupling coefficient $g_{\mathrm{2D}} = g \times (\int |\psi_0|^2 \psi_0 dz)/(\int \psi_0
dz)$.

%We have found by manipulating equations (\ref{variat1}) and (\ref{variat2}), that collapse of the
%soliton, when it occurs (at time $t_{col}$), is such that $W \approx V \sim (t_{col}-t)^{0.4}$, or
%else collapse is only in $\rho$ and $W \sim (t_{col}-t)^{0.5}$, $V \rightarrow V_{col}$ and
%$V-V_{col} \sim (t_{col}-t)\ln(t_{col}-t)$. In the 2D model, in which $V=V_0$, again $W \sim
%(t_{col}-t)^{0.5}$. This is an indication of one of the two above 3D scenarios. On the other hand,
%spreadout, when it occurs, is such that $W$ and/or $V$ are proportional to $t$ at large times.
%Alternatively, just $W \sim t$ and $V_0 \rightarrow$ const. A more detailed analysis will be
%presented in a fuller version of this work \cite{future}.

We simulated both the full GPE, Eq.~(\ref{NLS}), using an axisymmetric code (for 3D), a
Cartesian code (for 2D), and the variational equations for comparison. Numerical
simulations followed the path outlined in Fig.~\ref{fig1}. The parameters used in the
simulations would correspond, for $^{85}$Rb atoms and an OL period of $\lambda = 3\,\mu
$m, to an initial radial-confinement frequency of $\omega _{\perp }=2 \pi \times 39$ Hz,
a FRM frequency of $\Omega=2 \pi \times 2.87$ kHz. A good candidate is the Feshbach
resonance at 155 G for $^{85}$Rb. Since in this case the $a=0$ point is far from the
resonance, both two and three body losses are negligible \cite{155G}. A lattice depth of
$\varepsilon_{\mathrm{f}} =10.25$ recoil energies, and an effective nonlinear coefficient
of $(N/n)a=\pm \, 10^{-5}\,$m are required. Here $N/n$ is the number of atoms per lattice
cell, with a {\it total} number of atoms in the range of $10^{4}-10^{6}$. The respective
values of the normalized parameters are given in the figure captions.

\begin{figure}[tbp]
\includegraphics[width=8.5cm]{Dynamic2.eps}
\caption{Evolution of the amplitude of the central peak in the pattern comprising five cells of
the optical lattice (shown in the inset). The normalized parameters are $g_{0 \mathrm{f}}=-30$,
$g_{1 \mathrm{f}}=4g_{0 \mathrm{f}}$, $\varepsilon_{\mathrm{f}}=20.5$, $\Omega =22$, $\protect\omega _{\perp }=0.3$,
$t_{1}=30$, $t_{2}=100$, $t_{3}=120$, and $t_{4}=130$. The snapshots in the inset is taken at
$t=6000$. The unit of time is $m\lambda^2/(\pi^2\hbar$). \label{fig2}}
\end{figure}

Examples of the numerical results, which are generic, are displayed in Fig.~\ref{fig2},
which shows the evolution of the central-peak's amplitude versus time in dimensionless
units, defined by Eq.~(\ref{NLS}). After an initial transient, a stable structure is
established, featuring breathing without any systematic decay. The figure reveals the
influence of neighboring solitons on the amplitude of the central peak. The thin curve
corresponds to the evolution of the full multipeak structure, shown in the inset for a
fixed moment of time. To obtain the thick curve we repeated the above calculations up to
the time $t=7500$ and then removed all but the central peak. In this case the amplitude
gradually decreased. The interaction between neighboring solitons can be explained by the
difference between oscillation periods from cell to cell due to small deviations in the
numbers of atoms. It is essentially the same mechanism as in a Josephson junction. When
the lattice potential is weak, interaction between neighboring solitons is possible.
%The difference between these two cases proves that
%the dynamics of stable breathers as described here is ultimately a collective multi-peak
%phenomenon.

\begin{figure}[tbp]
\includegraphics[width=8.5cm]{ppsd233.eps}
\caption{(Color online). Stability regions for the 3D solitons in the $\left( |g_{0 \mathrm{f}}|,\Omega \right)$
plane, as predicted by the variational approximation (shaded area),
and found from direct simulations of the Gross-Pitaevskii equation (circles).
%Note that $g_{0 \mathrm{f}}$ is always negative.
Other parameters are as in Fig.~\protect\ref{fig2}, except $\varepsilon_{\mathrm{f}}=50$.
Note the similarity of the lower region to that of Fig.~\protect\ref{fig3b}. This region
corresponds to Q2D solitons. On the other hand, the upper region contains fully 3D
solitons. The borders of the VA stability regions were found analytically and will be
given in a fuller version \cite{future}. The lowest excitation frequency of the confining
potential is $\Omega_0=26.76$.\label{fig3a}}
\end{figure}

\begin{figure}[tbp]
\includegraphics[width=8.5cm]{ppsd232.eps}
\caption{(Color online). Same as Fig.~\protect\ref{fig3a}, but for Q2D treatment. \label{fig3b}}
\end{figure}

%In the 3D treatment we used both Eqs. (\ref{variat1}) and (\ref{variat2}).
In Figs.~\ref{fig3a} and \ref{fig3b} we have collected results of a systematic scan of
parameter space based on GPE simulations and compared them with the predictions of the VA
(a similar analysis can be performed if we replace 1D OL with a 1D harmonic trap - the
conclusion does not depend on the form of confinement in the $z$ direction). By stability
we mean shape preservation during one run. When this was not the case, after a short
period of time we observed clear collapse or spreadout in transverse directions.

The agreement between VA and direct simulations is very good. As seen from
Fig.~\ref{fig3a}, in our fully 3D treatment we found two islands of stability. Note the
similarity of the lower region to that of Fig.~\protect\ref{fig3b}, which portrays the
results of the Q2D treatment. This region corresponds to Q2D solitons. On the other hand,
the upper region contains fully 3D solitons. It appears when the frequency of modulation
exceeds the lowest excitation frequency of the confining potential, $\Omega_0=26.76$. If
the strength of the lattice $\varepsilon_{\mathrm{f}}$ is increased, this region moves
towards higher frequencies, the Q2D region expands upwards, and becomes more and more
like in Fig.~\protect\ref{fig3b}. This will be demonstrated in a fuller version of this
work \cite{future}.

%We compared the above results with the 2D treatment. When there is no modulation  the 2D treatment
%can predict the result of the 3D stability analysis, provided that the confinement is tight
%enough. However here the stability regions for the 2D and 3D problems are different .
%For moderate confinement the 3D region of stability is {\it very different} from that of
%Q2D (see Fig.~\ref{fig3a} and \ref{fig3b}). A large,
As we saw, new region of stability appears in the 3D treatment as compared to the Q2D
treatment. We are more used to the effect of adding a new dimension simply shrinking or
abolishing the basin of stability of solitons or waves. For example, water waves are
unstable with respect to perturbations along their direction of propagation only when the
depth exceeds a critical value \cite{TB}. When, however, two dimensional perturbations
are allowed, there will always be an unstable angle regardless of the depth \cite{WD}.
Another example is that of 1D solitons of the Nonlinear Schr\"{o}dinger equation (NLS)
with constant coefficients. They are stable in 1D, but unstable in 2D or 3D. This is also
true for some NLS waves \cite{EI}.

In the problem treated in this Letter this is not the case. For much of the quasi two
dimensional stability chart adding a degree of freedom {\it stabilizes} the soliton
solution. The key to this dichotomy seems to be the presence of a periodic modulation,
absent in the above mentioned classical examples. This can be illustrated by a simple
case involving oscillators. Take as the one dimensional version a forced oscillator
problem:
\begin{equation}
\ddot{x}+\omega_0^2x=y\cos(\omega_0t).
\end{equation}
If $y$ is fixed, the solution has a secular component $x=yt/(2\omega_0)\sin(\omega_0t) +
F(t)$, where $F(t)$ is a periodic function, and so the amplitude will grow as $t$. If
however we allow a second degree of freedom, such that $y$ also oscillates (for instance
$\ddot{y}+\epsilon^2y=0$) the solution stabilizes, unless $\epsilon = \pm 2\omega_0$. In
general this can be the case when there are periodic modulations. This fact, obvious in
oscillator theory, is perhaps less well known in the soliton context.

%Thus, adding a degree of freedom can stabilize, as in our main physical problem.
%In this work we compare regions of 2D and 3D stability of Feshbach Resonance Managed solitons in
%BECs where the interatomic interaction is modulated and on average attractive. We find that the 3D
%treatment has a region in parameter space which is not a subset of that of the quasi 2D treatment;
%we found regions of 3D stability, corresponding to instability in quasi 2D treatment. Therefore,
%we concluded that a quasi 2D treatment, convenient as it is  in predicting a general behavior of
%the system in the strong confinement limit, cannot provide a precise indication for real 3D
%stability (it is more successful in predicting the collapse phenomenon). This result is
%illustrated by a simple analogy.
%In previous work, 3D stability was inferred from the 2D treatment on the basis of just one example
%\cite{6Saito,9VPG}. The whole regions of stability must be considered before drawing any general conclusions.
%The possibility of reaching the result was demonstrated within the framework of the variational
%approximation (VA) and verified in direct simulations. The boundaries of the stability region
%obtained from the VA are confirmed by analytical considerations \cite{future}. The VA predictions
%are also strongly corroborated by the simulations. Thus, the VA proves to be very efficient in
%predicting where to look for the stability regions.

The main result of this letter is the possibility of creating fully 3D breather solitons
in a BEC confined by a 1D optical lattice potential, corresponding to the upper region in
Fig.~\ref{fig3a}. The stable patterns may feature a multi-cell structure, which in the
case studied here forms a set of weakly interacting fundamental solitons. The scheme
proposed in this Letter is based on a combination of FRM and a 1D optical lattice, and
could be implemented in an experiment, as outlined here. This would open the way to the
creation of robust 3D solitons (breathers) in BECs. A similar idea can apply in the field
of nonlinear optics. The Feshbach resonance could be replaced by a nonlinear periodic
structure, for example ferromagnetic domains as used in contemporary quasi phase matching
\cite{Arie}, or by optically induced photorefractive lattices \cite{Wiesio}.

M.M. acknowledges support from the KBN grant 2P03 B4325, M.T. was supported by the Polish Ministry
of Scientific Research and Information Technology under grant PBZ MIN-008/P03/2003 and E.I. from
that of grant 2P03B09722. The work of B.A.M. was partially supported by the Israel Science
Foundation through grant No. 8006/03.

\begin{thebibliography}{99}

\bibitem{1fsol} S.~Burger \textit{et al.},
%%K.~Bongs, S.~Dettmer, W.~Ertmer, and K.~Sengstock,
Phys. Rev. Lett. \textbf{83}, 5198 (1999); J.~Denschlag \textit{et al.},
%%J.~E.~Simsarian, D.~L.~Feder, C. W. Clark, L. A. Collins,
%%J. Cubizolles, L. Deng, E. W. Hagley, K. Helmerson, W. P. Reinhardt,
%%S. L. Rolston, B. I. Schneider, and W. D. Phillips,
Science \textbf{287}, 97 (2000).

\bibitem{2bright} L.~Khaykovich \textit{et al.},
%%F. Scherck, G. Ferrari, T. Bourdel,
%%J.Cubizolles, L. D. Carr, Y. Castin, C. Salomon,
Science \textbf{296}, 1290 (2002), K.~E.~Strecker \textit{et al.},
%%G. B. Partridge, A. G. Truscott and R. G. Hulet,
Nature \textbf{417}, 153 (2002).

\bibitem{ober} B.~Eiermann \textit{et al.},
%Th. Anker, M. Albiez, M. Taglieber, P. Treutlein, K.-P. Marzlin, and M. K. Oberthaler,
Phys. Rev. Lett. \textbf{92}, 230401 (2004).

\bibitem{3Inouye} S.~Inouye \textit{et al.},
%%M. R. Andrews, J. Stenger, H.-J. Miesner,
%%D. M. Stamper-Kurn, and W. Ketterle,
Nature \textbf{392}, 151 (1998); E.~A.~Donley \textit{et al.},
%%N. R. Claussen, S. L. Cornish, J. L. Roberts, E. A. Cornell,
%%and C. E. Wieman,
Nature \textbf{412}, 295 (2001); H.~Saito and M.~Ueda, Phys. Rev. A \textbf{65}, 033624 (2002).
%;P. O.~Fedichev \textit{et al.},
%Yu. Kagan, G. V. Shlyapnikov, and J. T. M. Walraven,
%Phys. Rev. Lett. \textbf{77}, 2913 (1996); M.~Theis \textit{et al.},
%G. Thalhammer, K. Winkler, M. Hellwig, G. Ruff, R. Grimm, and J.
%H. Denschlag,
%Phys. Rev. Lett. \textbf{93}, 123001 (2004),

\bibitem{4Greece} P.~G.~Kevrekidis \textit{et al.}
%G. Theocharis, D. J. Frantzeskakis, and B. A. Malomed,
Phys. Rev. Lett. \textbf{90}, 230401 (2003).

%\bibitem{5Isaac} I. Towers and B.A. Malomed, J. Opt. Soc. Am. \textbf{19},
%537 (2002).

%\bibitem{7Carr} L. D. Carr and Y. Castin, Phys. Rev. A \textbf{66}, 063602
%(2002).

\bibitem{6Saito} H.~Saito and M.~Ueda, Phys. Rev. Lett. \textbf{90}, 040403
(2003).

\bibitem{8Abdul} F.~Abdullaev \textit{et al.},
%J. G. Caputo, R. A. Kraenkel, and B. A. Malomed,
Phys. Rev. \textbf{A 67}, 013605 (2003).

\bibitem{9VPG} G.~D.~Montesinos, V.~M.~Perez-Garcia, and P.~J.~Torres,
Physica D \textbf{191}, 193 (2004).

\bibitem{Gerlitz} A.~G\"orlitz {\it et al},
%J.M.Vogels, A. E. Leanhardt, C. Raman, T. L.
%Gustavson, J. R. Abo-Shaeer, A. P. Chikkatur, S. Gupta, S. Inouye, T. Rosenband,
%and W. Ketterle,
Phys. Rev. Lett. \textbf{87}, 130402 (2001).

\bibitem{10EL} M.~Trippenbach, M.~Matuszewski and B.~A.~Malomed,
Europhys. Lett. (in press), quant-ph/0410235.

\bibitem{11OL} H.~Stecher \textit{et al.},
%%H. Ritsch, P. Zoller, F. Sander, T. Esslinger, and T. W. Hansch,
Phys. Rev. A \textbf{55}, 545 (1997).

\bibitem{12BBB} B.~B.~Baizakov, B.~A.~Malomed, and M.~Salerno, Europhys. Lett.
\textbf{63}, 642 (2003).

\bibitem{Progress} B.~A.~Malomed, Progr. Opt. \textbf{43}, 69 (2002).

\bibitem{155G} J.~L.~Roberts, N.~R.~Claussen, S.~L.~Cornish, and C.~E.~Wieman
Phys. Rev. Lett. \textbf{85}, 728 (2000).

\bibitem{future}  M.~Matuszewski {\it et al.} (in preparation).

\bibitem{TB} T.~B.~Benjamin and J.~E.~Feir, J. Fluid. Mech. \textbf{27}, 417 (1967).

\bibitem{WD} W.~D.~Hayes, Proc. R. Soc. London, \textbf{A 332}, 199 (1973).

\bibitem{EI} E.~Infeld and G.~Rowlands, \textit{Nonlinear Waves, Solitons and Chaos}
(Cambridge University Press, Cambridge, 2nd edition, 2000).
%; D.~Anderson, A.~Bondeson and M.~Lisak,
%J. Plasma Phys. \textbf{21}, 259 (1979); E.~Infeld and G.~Rowlands, Z. Phys. \textbf{B 37}, 277
%(1980).

%% \bibitem{proza} To create stable 3D solitons in an experiment, one needs an OL of
%% strength $\varepsilon_{\mathrm{f}}$ of at least $10
%% E_{\mathrm{recoil}}$. We define $\varepsilon_{\mathrm{rel}}=2 \varepsilon_{\mathrm{f}} /
%% E_{\mathrm{recoil}}$. The modulation frequency $\Omega$ must be in the range
%% $[\Omega_0,2\Omega_0]$, where $\Omega_0=\frac{2\pi^2 \hbar}{m \lambda^2}
%% \sqrt{4\varepsilon_{\mathrm{rel}}- \sqrt{\varepsilon_{\mathrm{rel}}}}$ is the excitation frequency
%% of the confining potential. The final value of the mean scattering length $\overline{a_0}$ and
%% number of atoms $N$ must satisfy the condition  $\overline{a_0}<a_c=-n
%% \lambda\varepsilon_{\mathrm{rel}}^{1/4}/(2 \sqrt{\pi}N)$. The scattering length modulations $a_1$
%% should be a few times larger than $|\overline{a_0}|$. Finally the initial radial optical trap
%% frequency $\omega_{\perp}$ should be smaller than $\omega_{\perp 0}=m \lambda^3
%% (\frac{\overline{a_0}}{a_c}-1)(\Omega^2-\Omega_0^2)/(2 N a_1 \varepsilon_{\mathrm{rel}} \hbar
%% \pi^{5/2})$.

\bibitem{Arie} G. Rosenman et al, Appl. Phys. Lett. 73, 3650 (1998).

\bibitem{Wiesio} D. Neshev et al, Phys. Rev. Lett. 93, 083905 (2004).

\end{thebibliography}

\end{document}

cond-mat0506545.extracted_bib
\begin{thebibliography}{99}

\bibitem{1fsol} S.~Burger \textit{et al.},
%%K.~Bongs, S.~Dettmer, W.~Ertmer, and K.~Sengstock,
Phys. Rev. Lett. \textbf{83}, 5198 (1999); J.~Denschlag \textit{et al.},
%%J.~E.~Simsarian, D.~L.~Feder, C. W. Clark, L. A. Collins,
%%J. Cubizolles, L. Deng, E. W. Hagley, K. Helmerson, W. P. Reinhardt,
%%S. L. Rolston, B. I. Schneider, and W. D. Phillips,
Science \textbf{287}, 97 (2000).

\bibitem{2bright} L.~Khaykovich \textit{et al.},
%%F. Scherck, G. Ferrari, T. Bourdel,
%%J.Cubizolles, L. D. Carr, Y. Castin, C. Salomon,
Science \textbf{296}, 1290 (2002), K.~E.~Strecker \textit{et al.},
%%G. B. Partridge, A. G. Truscott and R. G. Hulet,
Nature \textbf{417}, 153 (2002).

\bibitem{ober} B.~Eiermann \textit{et al.},
%Th. Anker, M. Albiez, M. Taglieber, P. Treutlein, K.-P. Marzlin, and M. K. Oberthaler,
Phys. Rev. Lett. \textbf{92}, 230401 (2004).

\bibitem{3Inouye} S.~Inouye \textit{et al.},
%%M. R. Andrews, J. Stenger, H.-J. Miesner,
%%D. M. Stamper-Kurn, and W. Ketterle,
Nature \textbf{392}, 151 (1998); E.~A.~Donley \textit{et al.},
%%N. R. Claussen, S. L. Cornish, J. L. Roberts, E. A. Cornell,
%%and C. E. Wieman,
Nature \textbf{412}, 295 (2001); H.~Saito and M.~Ueda, Phys. Rev. A \textbf{65}, 033624 (2002).
%;P. O.~Fedichev \textit{et al.},
%Yu. Kagan, G. V. Shlyapnikov, and J. T. M. Walraven,
%Phys. Rev. Lett. \textbf{77}, 2913 (1996); M.~Theis \textit{et al.},
%G. Thalhammer, K. Winkler, M. Hellwig, G. Ruff, R. Grimm, and J.
%H. Denschlag,
%Phys. Rev. Lett. \textbf{93}, 123001 (2004),

\bibitem{4Greece} P.~G.~Kevrekidis \textit{et al.}
%G. Theocharis, D. J. Frantzeskakis, and B. A. Malomed,
Phys. Rev. Lett. \textbf{90}, 230401 (2003).

%\bibitem{5Isaac} I. Towers and B.A. Malomed, J. Opt. Soc. Am. \textbf{19},
%537 (2002).

%\bibitem{7Carr} L. D. Carr and Y. Castin, Phys. Rev. A \textbf{66}, 063602
%(2002).

\bibitem{6Saito} H.~Saito and M.~Ueda, Phys. Rev. Lett. \textbf{90}, 040403
(2003).

\bibitem{8Abdul} F.~Abdullaev \textit{et al.},
%J. G. Caputo, R. A. Kraenkel, and B. A. Malomed,
Phys. Rev. \textbf{A 67}, 013605 (2003).

\bibitem{9VPG} G.~D.~Montesinos, V.~M.~Perez-Garcia, and P.~J.~Torres,
Physica D \textbf{191}, 193 (2004).

\bibitem{Gerlitz} A.~G\"orlitz {\it et al},
%J.M.Vogels, A. E. Leanhardt, C. Raman, T. L.
%Gustavson, J. R. Abo-Shaeer, A. P. Chikkatur, S. Gupta, S. Inouye, T. Rosenband,
%and W. Ketterle,
Phys. Rev. Lett. \textbf{87}, 130402 (2001).

\bibitem{10EL} M.~Trippenbach, M.~Matuszewski and B.~A.~Malomed,
Europhys. Lett. (in press), quant-ph/0410235.

\bibitem{11OL} H.~Stecher \textit{et al.},
%%H. Ritsch, P. Zoller, F. Sander, T. Esslinger, and T. W. Hansch,
Phys. Rev. A \textbf{55}, 545 (1997).

\bibitem{12BBB} B.~B.~Baizakov, B.~A.~Malomed, and M.~Salerno, Europhys. Lett.
\textbf{63}, 642 (2003).

\bibitem{Progress} B.~A.~Malomed, Progr. Opt. \textbf{43}, 69 (2002).

\bibitem{155G} J.~L.~Roberts, N.~R.~Claussen, S.~L.~Cornish, and C.~E.~Wieman
Phys. Rev. Lett. \textbf{85}, 728 (2000).

\bibitem{future}  M.~Matuszewski {\it et al.} (in preparation).

\bibitem{TB} T.~B.~Benjamin and J.~E.~Feir, J. Fluid. Mech. \textbf{27}, 417 (1967).

\bibitem{WD} W.~D.~Hayes, Proc. R. Soc. London, \textbf{A 332}, 199 (1973).

\bibitem{EI} E.~Infeld and G.~Rowlands, \textit{Nonlinear Waves, Solitons and Chaos}
(Cambridge University Press, Cambridge, 2nd edition, 2000).
%; D.~Anderson, A.~Bondeson and M.~Lisak,
%J. Plasma Phys. \textbf{21}, 259 (1979); E.~Infeld and G.~Rowlands, Z. Phys. \textbf{B 37}, 277
%(1980).

%% \bibitem{proza} To create stable 3D solitons in an experiment, one needs an OL of
%% strength $\varepsilon_{\mathrm{f}}$ of at least $10
%% E_{\mathrm{recoil}}$. We define $\varepsilon_{\mathrm{rel}}=2 \varepsilon_{\mathrm{f}} /
%% E_{\mathrm{recoil}}$. The modulation frequency $\Omega$ must be in the range
%% $[\Omega_0,2\Omega_0]$, where $\Omega_0=\frac{2\pi^2 \hbar}{m \lambda^2}
%% \sqrt{4\varepsilon_{\mathrm{rel}}- \sqrt{\varepsilon_{\mathrm{rel}}}}$ is the excitation frequency
%% of the confining potential. The final value of the mean scattering length $\overline{a_0}$ and
%% number of atoms $N$ must satisfy the condition  $\overline{a_0}<a_c=-n
%% \lambda\varepsilon_{\mathrm{rel}}^{1/4}/(2 \sqrt{\pi}N)$. The scattering length modulations $a_1$
%% should be a few times larger than $|\overline{a_0}|$. Finally the initial radial optical trap
%% frequency $\omega_{\perp}$ should be smaller than $\omega_{\perp 0}=m \lambda^3
%% (\frac{\overline{a_0}}{a_c}-1)(\Omega^2-\Omega_0^2)/(2 N a_1 \varepsilon_{\mathrm{rel}} \hbar
%% \pi^{5/2})$.

\bibitem{Arie} G. Rosenman et al, Appl. Phys. Lett. 73, 3650 (1998).

\bibitem{Wiesio} D. Neshev et al, Phys. Rev. Lett. 93, 083905 (2004).

\end{thebibliography}