Dynamics of Bose–Einstein condensates near Feshbach resonance in external potential

Abstract

We review our recent theoretical advances in the dynamics of Bose–Einstein condensates with tunable interactions using Feshbach resonance and external potential. A set of analytic and numerical methods for Gross–Pitaevskii equations are developed to study the nonlinear dynamics of Bose–Einstein condensates. Analytically, we present the integrable conditions for the Gross–Pitaevskii equations with tunable interactions and external potential, and obtain a family of exact analytical solutions for one- and two-component Bose–Einstein condensates in one and two-dimensional cases. Then we apply these models to investigate the dynamics of solitons and collisions between two solitons. Numerically, the stability of the analytic exact solutions are checked and the phenomena, such as the dynamics and modulation of the ring dark soliton and vector-soliton, soliton conversion via Feshbach resonance, quantized soliton and vortex in quasi-two-dimensional are also investigated. Both the exact and numerical solutions show that the dynamics of Bose–Einstein condensates can be effectively controlled by the Feshbach resonance and external potential, which offer a good opportunity for manipulation of atomic matter waves and nonlinear excitations in Bose–Einstein condensates.

Graphical abstract

Keywords

Bose–Einstein condensate / Feshbach resonance / soliton

1 Introduction

Solitons are fundamental excitations of nonlinear media and have attracted great interest in diverse contexts of science and engineering, such as the dynamics of waves in shallow water, transport along DNA and other macromolecules, plasmas, fiber optic communications and cold atoms. The realization of Bose–Einstein condensates (BECs) in gases of weakly interacting alkali-metal atoms [1, 2] introduces an unparalleled platform for the study of these nonlinear excitations, where both bright and dark solitons have been observed [3-6] and the gap solitons are realized in an optical lattice experimentally [7, 8]. These offer an unusual opportunity to study the dynamics of BECs and solitons, and the relevant topics.

Moreover, the important mechanism, Feshbach resonance, provides a unique tool to further control the BECs system. It is well known that interatomic interactions greatly affect a number of properties of BECs, including both static and dynamic ones [9-12]. The Feshbach resonance makes it possible to change the interaction strength, even its sign, with a tunable time dependent magnetic, or optical field. This degree of freedom offers us a good opportunity for the manipulation of atomic matter waves and nonlinear excitations in BECs [13-20]. On the theoretical side, several forms of the time-dependent scattering lengths, such as the exponential function [21-23], and the periodic function [24-28], have been proposed and studied. Here, we want to point out that the condensates with a spatially modulated nonlinearity by manipulating scattering length locally have been proposed. In this case, the spatial dependence of scattering length can be implemented by a spatially inhomogeneous external magnetic field in the vicinity of a Feshbach resonance [29].

With the development of experiment, two (or more) internal states or different atoms can be populated. This degree of freedom has stimulated considerable attention in vector solitons [30-36], and has given rise to rich phenomena and complex dynamics in the two-component BECs such as soliton trains, soliton pairs, multidomain walls, and multimode collective [37-43]. Vector solitons consist of more than one field components, which are supported not only by self-interaction of the same fields, but also by cross interactions between fields of different type. Recently, experimental observation of heteronuclear Feshbach molecules from a ⁸⁷Rb and ⁸⁵Rb gas, and tunable miscibility in a dual-component Bose–Einstein condensates have been realization in [44]. Meanwhile, the vector solitons in the form of “dark–dark” [45], “bright–bright” solitons [46] and “bright–dark” [47]are studied.

In this article, we consider the dynamics of BECs confined in an external trap with longitudinal and transverse confining frequencies (denoted by {\omega }_z and {\omega }_{\perp }, respectively). If {\omega }_z\ll {\omega }_{\perp }, the system can be regarded as in the quasi-one-dimensional (1D) geometry; On the other hand, the system can be regarded as in the 2D geometry if {\omega }_{\perp }\ll {\omega }_z. The scope of this brief review does not allow us to cover all the important issues above. Instead, we focus our attention on a few topics about the dynamics of BECs near Feshbach resonance in external potential. In Section 2, we consider the dynamics of quasi-1D BECs with feeding parameters or not. Section 3 is devoted to the discussion of dynamics of 2D BECs. Finally, we summarize this review.

2 Dynamics of BECs in one-dimension

2.1 One-component BECs in an expulsive parabolic potential

At the mean-field level, condensate dynamics at zero temperature obey a nonlinear wave equation known as the Gross–Pitaevskii equation

\begin{array}{}\text{(1)}& \mathrm{i}\mathrm{\hslash }{\displaystyle \frac{\partial \psi }{\partial t}}=-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2m}}{\nabla }^2\psi +V_{\mathrm{ext}}\left(\mathit{r}\right)\psi +g{\left|\psi \right|}^2\psi \end{array}

where g=\frac{4\pi {\mathrm{\hslash }}^2{a}_s}{m}with a_s the s-wave scattering length. In the case of cigar-shaped BECs, where {\omega }_y={\omega }_z={\omega }_{\perp } and \lambda =2{\omega }_x/{\omega }_{\perp }\ll 1, it is reasonable to reduce the GP equation into an effective 1D dimensionless nonlinear Schrödinger equation

\begin{array}{}\text{(2)}& \mathrm{i}\frac{\partial \psi (x,t)}{\partial t}+\frac{{\partial }^2\psi (x,t)}{\partial x^2}+2a(t){|\psi (x,t)|}^2\psi (x,t)+\frac{1}{4}{\lambda }^2{x}^2\psi (x,t)=0\end{array}

where time t and coordinate x are measured in units 2/{\omega }_{\perp } and a_{\perp }, where a_{\perp }={(\mathrm{\hslash }/m{\omega }_{\perp })}^{1/2} and a_0={(\mathrm{\hslash }/m{\omega }_1)}^{1/2} are linear oscillator lengths in the transverse and cigar-axis directions, respectively. {\omega }_{\perp } and {\omega }_1 are respective harmonic oscillator frequencies, and \lambda =2|{\omega }_1|/{\omega }_{\perp }\ll 1.

When the scattering length is increased in the form of a(t)=g_0\exp (\lambda t), we perform the Darboux transformation to obtain the exact solution of Eq. (2) as follows:

\begin{array}{}\text{(3)}& \psi =[A_0+A_1\frac{(\gamma \cosh \theta +\cos \phi )+\mathrm{i}(\alpha \sinh \theta +\beta \sin \phi )}{\cosh \theta +\gamma \cos \phi }]\cdot \exp (\frac{\lambda t}{2}+\mathrm{i}{\phi }_c)\end{array}

where

\begin{array}{}\text{(3)}& \begin{array}{c}\theta =-{\displaystyle \frac{\left[\left(k_0+k_s\right){\Delta }_R-\sqrt{g_0}{A}_1{\Delta }_I\right]\left[\exp \left(2\lambda t\right)-1\right]}{2\lambda }}+{\Delta }_{R}x\cdot \exp \left(\lambda t\right)\hfill \\ \phi =-{\displaystyle \frac{\left[\left(k_0+k_s\right){\Delta }_I+\sqrt{g_0}{A}_1{\Delta }_R\right]\left[\exp \left(2\lambda t\right)-1\right]}{2\lambda }}+{\Delta }_{I}x\exp \left(\lambda t\right)\hfill \\ \alpha ={\displaystyle \frac{\sqrt{g_0}{A}_0\left(k_0-k_s+{\Delta }_I\right)}{\Lambda }}\hfill \\ \beta =1-{\displaystyle \frac{2g_0{A}_0^2}{\Lambda }},\gamma ={\displaystyle \frac{\sqrt{g_0}{A}_0\left({\Delta }_R-\sqrt{g_0}{A}_1\right)}{\Lambda }}\hfill \\ \Delta =\sqrt{{\left[-\sqrt{g_0}{A}_1+\mathrm{i}\left(k_s-k_0\right)\right]}^2-4g_0{A}_0^2}\equiv {\Delta }_R+\mathrm{i}{\Delta }_I\hfill \\ \Lambda =g_0{A}_0^2+{\displaystyle \frac{\left({\Delta }_R-\sqrt{g_0}{A}_1^2\right)}{4}}+{\displaystyle \frac{{\left(k_s-k_0+{\Delta }_I\right)}^2}{4}}\hfill \end{array}\end{array}

k_s is the arbitrary real constant. On the one hand, when A_0=k_0=0, Eq. (3) reduces to the well-known one solitonic solution. On the other hand, when the amplitude of this soliton vanishes (A_1=0), Eq. (3) reduces to plane wave solution. Thus, Eq. (3) represents a bright soliton embedded in the background of density.

In the following discussion, we found that the manipulation of the scattering length can be used to compress a bright soliton of BECs into an assumed peak matter density. For simplicity, we assume k_0=k_s and only consider the case of A_1^2>4A_0^2. In this case, Eq. (3) can be deduced as:

\begin{array}{}\text{(4)}& \psi =(-A_0+{\delta }_2\frac{{\delta }_2\cos \phi -\mathrm{i}{A}_1\sin \theta }{A_1\cosh \theta -2A_0\cos \phi })\cdot \exp (\frac{\lambda t}{2}+\mathrm{i}{\phi }_c)\end{array}

where \phi =-\frac{\sqrt{g_0}{A}_1{\delta }_2\left[\exp (2\lambda t)-1\right]}{2\lambda },{\delta }_2=\sqrt{A_1^2-4A_0^2}, and \theta =\sqrt{g_0}{\delta }_{2}x\exp (\lambda t)-\frac{\sqrt{g_0}{k}_0{\delta }_2\left[\exp (2\lambda t)-1\right]}{\lambda }.

Figure 1 shows the density distribution of BECs in the expulsive parabolic potential. With time going on, the bright soliton has an increase in the peaking value and a compression in the width. With this advantage, we can get a bright soliton with the assumed peak matter density. This degree of freedom provides us an experimental tool for investigating the range of validity of the 1D GP equation.

Furthermore, we find that when \sinh \theta =0, the peak matter density of bright soliton arrives at the maximum; while when \cosh \theta =\frac{A_1}{A_0\cos \phi }-\frac{A_0\cos \phi }{A_1}, the peak matter density of bright soliton arrives at the minimum. This means that the bright soliton can only be squeezed into the assumed peak matter density between the minimum and maximum values. In order to investigate the stability of the bright soliton, we obtain {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\left[{\left|\psi (x,t)\right|}^2-{\left|\psi (\pm \mathrm{\infty },t)\right|}^2\right]dx=\frac{2{\delta }_2}{\sqrt{g_0}}, which is the exact number of the atoms in the bright soliton. Thus during the process of the compression of the bright soliton, the number of atoms in the bright soliton stays invariant. On the other hand, the quantity \kappa ={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\left|\psi (x,t)-\psi (\pm \mathrm{\infty },t)\right|}^{2}dx=\frac{2{\delta }_2}{\sqrt{g_0}}\left(1+A_{0}M\cos \phi \right) with M=4\arctan \frac{\sqrt{A_1+2A_0\cos \phi }}{\sqrt{A_1-2A_0\cos \phi }}/\sqrt{A_1^2-4A_0^2{\cos }^2\phi } displays that a time-periodic atomic exchange is formed between the bright soliton and the background. As shown in Fig. 2(a), in the case of zero-background, there will not be an exchange of atoms. While in the case of nonzero-background [Fig. 2(b) and (c)], the exchange of atoms between the bright soliton and the background becomes more quick with the increase of the absolute value of the scattering length [23].

2.2 The integrable condition

We also start by considering the quasi-1D nonlinear Schrödinger equation

\begin{array}{}\text{(5)}& \mathrm{i}\frac{\partial \psi }{\partial t}=-\frac{1}{2}\frac{{\partial }^2\psi }{\partial x^2}+\frac{a_s(t)}{a_{\mathrm{B}}}{|\psi |}^2\psi +\frac{{\omega }_1^2}{2{\omega }_{\perp }^2}{x}^2\psi \end{array}

If we allow the axial frequency of the harmonic trap to become also time dependent {\omega }_1={\omega }_1(t), and require it to satisfy the following integrability relation with the scattering length a_s(t):

\begin{array}{}\text{(6)}& -{\displaystyle \frac{1}{a_s(t)}}\frac{{\mathrm{d}}^2{a}_s(t)}{\mathrm{d}{t}^2}+\frac{2}{a_s^2(t)}{\left[\frac{\mathrm{d}{a}_s(t)}{\mathrm{d}t}\right]}^2+\frac{{\omega }_1^2(t)}{{\omega }_{\perp }^2}=0\end{array}

Then, the nonlinear GP equation (5) can be reduced to the standard nonlinear Schrödinger equation and exactly solved, with the following general solution:

\begin{array}{}\text{(7)}& \psi (x,t)=\exp [-{\displaystyle {\int }_{t_0}^t}\Gamma (t)dt]\varphi (X,T)\exp [\mathrm{i}\Gamma (t)x^2]\end{array}

where φ is the well-known solution for the standard Schrödinger equation, with new spatial and temporal variables X=A_1\exp [-2{\int }_{t_0}^t\Gamma (t)dt]x,T=\frac{A_1^2}{2}\cdot {\int }_{t_0}^t\exp [-4{\int }_{t_0}^t\Gamma (t^{\prime })d{t}^{\prime }]dt.A_1 is a real constant, which together with \Gamma (t) are determined by

\begin{array}{}\text{(8)}& a_s(t)=\pm a_{\mathrm{B}}{A}_1^2\exp \left[-2{\displaystyle {\int }_{t_0}^t}\Gamma (t)dt\right]\end{array}

From the above equations, we see that once the scattering length is determined, the trap frequency can also be determined. Furthermore, we find three special cases for the integrable conditions, namely, linear, exponential and sinusoidal time dependence of the magnetic field via Feshbach resonance. Note that the exact solution \psi (x,t) can be obtained for arbitrary time dependence of a_s(t), since we can always choose an appropriate time-dependent axial frequency, {\omega }_1(t), to satisfy the integrability relation [48, 49].

2.3 Dynamics of BECs in one dimension complex potential

In real experiments, atom loss and feeding of condensate from the thermal cloud are inevitable courses. In this subsection, we consider a one-dimension nonlinear Schrödinger equation with time-dependent atomic scattering length in an expulsive parabolic and complex potential

\begin{array}{}\text{(9)}& \mathrm{i}\frac{\partial \psi }{\partial t}=-\frac{{\partial }^2\psi }{\partial x^2}+2a(t){|\psi |}^2\psi -\frac{1}{4}{\lambda }^2{x}^2\psi +\mathrm{i}\gamma \psi \end{array}

where γ is a small parameter related to the feeding of condensate from the thermal cloud [50, 51]. In the following discussion, we develop a direct method to derive two families of exact solitons of Eq. (9), then give some thorough analysis for a bright soliton, a train of bright solitons and a dark soliton. Our results show that the lifetime of a bright or a dark soliton in BECs can keep longer times by reducing both the ratio of the axial oscillation frequency to radial oscillation frequency and the loss of atoms. It is demonstrated that a train of bright solitons in BECs may be excited with a strong enough background.

We first assume the solutions of Eq. (9) as follows:

\begin{array}{}\text{(10)}& \psi =[A_0(t)+A_1(t)\frac{\delta \cosh \xi +\cos \eta }{\cosh \xi +\delta \cos \eta }+\mathrm{i}{B}_1(t)\frac{\alpha \sinh \xi +\beta \sin \eta }{\cosh \xi +\delta \cos \eta }]\exp (\mathrm{i}\Delta )\end{array}

where \Delta =k_0(t)+k_1(t)x+k_2(t)x^2,\xi =p_1(t)x+q_1(t),\eta =p_2(t)x+q_2(t), and A_0(t),A_1(t),B_1(t),p_1(t),q_1(t),p_2(t),q_2(t),k_0(t),k_1(t),k_2(t) are real functions of t to be determined, and \alpha ,\beta ,\delta are real constants.

Substituting Eq. (10) into Eq. (9) we derive a set of ordinary differential equations (ODEs) with respect to a(t),A_0(t),A_1(t),B_1(t),p_1(t),q_1(t),p_2(t),q_2(t),k_0(t),k_1(t),k_2(t). Finally, solving these ODEs, we can obtain two families of analytical solutions of Eq. (9).

When interaction between atoms is attractive, the solution of Eq. (9) can be written as:

\begin{array}{}\text{(11)}& {\psi }_1=\Omega (A_c+A_s\frac{\delta \cosh \xi +\cos \eta }{\cosh \xi +\delta \cos \eta }+\mathrm{i}{A}_s\frac{\alpha \sinh \xi +\beta \sin \eta }{\cosh \xi +\delta \cos \eta })\exp (\mathrm{i}\Delta +\gamma t)\end{array}

where

\begin{array}{}\text{(11)}& \begin{array}{c}\Delta =k_2(t)x^2+k_1{\Omega }^{2}x+(2g_0{A}_c^2-k_1^2){\displaystyle \int {\Omega }^{4}dt}\hfill \\ \xi =\sqrt{g_0}{A}_s\{\beta {\Omega }^{2}x-2[k_1\beta +p_2-\frac{2g_0{p}_2{A}_c^2}{(g_0{A}_s^2+p_2^2)}]{\displaystyle \int }{\Omega }^4\mathrm{d}t\}\hfill \\ \eta =p_2{\Omega }^{2}x-[2p_2{k}_1+(p_2^2-g_0{A}_s^2)\beta ]{\displaystyle \int {\Omega }^{4}dt}\hfill \\ a(t)=-g_0{\Omega }^2\exp (-2\gamma t),\alpha =-\frac{2\sqrt{g_0}{A}_c{p}_2}{g_0{A}_s^2+p_2^2}\hfill \\ {\beta }^2=\frac{p_2^2+g_0(A_s^2-4A_c^2)}{g_0{A}_s^2+p_2^2},\delta =-\frac{2g_0{A}_c{A}_s}{g_0{A}_s^2+p_2^2}\hfill \\ \Omega =\exp [{\displaystyle \int }-2k_2(t)\mathrm{d}t],k_2(t)=\{\pm \frac{\lambda }{4},\frac{\lambda }{4}\tanh (\lambda t)\}\hfill \end{array}\end{array}

and A_c,A_s,g_0>0,p_2,k_1,\gammaare arbitrary real constants.

When interaction between atoms is repulsive, the solution of Eq. (9) can be written as:

\begin{array}{}\text{(12)}& {\psi }_2=\Omega (A_c+\mathrm{i}{A}_s\tanh \xi )\exp (\mathrm{i}\Delta +\gamma t)\end{array}

where

\begin{array}{}\text{(12)}& \begin{array}{c}\xi =\pm \sqrt{g_0}{\Omega }^{2}x+2A_s(\sqrt{g_0}{k}_1+g_0{A}_c){\displaystyle \int {\Omega }^{4}dt}\hfill \\ \Delta =k_2(t)x^2+k_1{\Omega }^{2}x-[2g_0(A_c^2+A_s^2)+k_1^2]{\displaystyle \int {\Omega }^{4}dt}\hfill \\ a(t)=g_0{\Omega }^2\exp (-2\gamma t)\hfill \end{array}\end{array}

and \Omega ,k_2(t) are the same as in (11), A_c,A_s,g_0>0,k_1 and γ are arbitrary real constants.

For the sake of simplicity, we will only discuss and analyze Eqs. (11) and (12) with k_2(t)=\lambda /[4\tanh (\lambda t)] and {\Omega }^2=\mathrm{sech}(\lambda t)/2.

Consider the dynamics of the bright soliton in the background. As shown in Fig. 3(a), when \alpha =p_2=0 under the realistic experiment parameters in Ref. [5], i.e., {\omega }_{\perp }=2\pi \times 710\mathrm{Hz},{\omega }_0=2\mathrm{i}\pi \times 70\mathrm{Hz},\lambda =-2|{\omega }_0|/{\omega }_{\perp }\approx -0.197,\gamma =-0.01 and g_0=0.4\mathrm{nm}, we can see that the lifetime of the BECs is about 20\times 4.5\times {10}^{-4}\mathrm{s}=9\mathrm{ms}. When \lambda =-0.02 (which can be derived from {\omega }_0=2\mathrm{i}\pi \times 7\mathrm{Hz} and {\omega }_{\perp }=2\pi \times 710\mathrm{Hz}) and \gamma =-0.001, from Fig. 3(b), the lifetime of the BECs can reach about 200 units of the dimensionless time, which corresponds to a real time of 0.1 s. Thus, in order to extending the lifetime of a soliton in BECs, we should take appropriate measures to reduce the absolute value of λ and γ.

When \beta =0, the solution {\psi }_1 is written as:

\begin{array}{}\text{(13)}& {\psi }_{12}=\Omega (A_c+A_s\frac{\delta \cosh \xi +\cos \eta +\mathrm{i}\alpha \sinh \xi }{\cosh \xi +\delta \cos \eta })\cdot \exp (\mathrm{i}\Delta +\gamma t)\end{array}

where 4A_c^2-A_s^2>0 and

\begin{array}{}\text{(13)}& \begin{array}{c}\xi =\sqrt{g_0}{A}_s{p}_2{\displaystyle \int {\Omega }^{4}dt},\alpha =-\frac{p_2}{2\sqrt{g_0}{A}_c}\hfill \\ \eta =p_2{\Omega }^{2}x-2p_2{k}_1{\displaystyle \int {\Omega }^{4}dt},\delta =-\frac{A_s}{2A_c}\hfill \\ p_2^2=g_0(4A_c^2-A_s^2),\Omega =\sqrt{\frac{\mathrm{sech}(\lambda t)}{2}}\hfill \end{array}\end{array}

It is easy to find that {\psi }_{12} is periodic with a period \Gamma =4\pi /[p_2\mathrm{sech}(\lambda t)] in the space coordinate and aperiodic in the temporal variable. As shown in Fig. 4, it is a train of bright solitons when \lambda =-0.197,A_c=7,A_s=12,g_0=0.4,k_1=\gamma =-0.01, thus we can conclude that an important condition for exciting a train of bright solitons is that the background is strong enough [50].

2.4 Dynamics of two-component BECs with time-dependent interactions

The mean-field dynamics of a two-species BECs is governed by the following equations:

\begin{array}{}\text{(14)}& \begin{array}{c}\mathrm{i}\frac{\partial {\psi }_1}{\partial t}=(-\frac{1}{2}\frac{{\partial }^2}{\partial x^2}+\frac{{\lambda }_1^2}{2}{x}^2+b_{11}{|{\psi }_1|}^2+b_{12}{|{\psi }_2|}^2){\psi }_1\hfill \\ \mathrm{i}\frac{\partial {\psi }_2}{\partial t}=(-\frac{\kappa }{2}\frac{{\partial }^2}{\partial x^2}+\frac{{\lambda }_2^2}{2\kappa }{x}^2+b_{21}{|{\psi }_1|}^2+b_{22}{|{\psi }_2|}^2){\psi }_2\hfill \end{array}\end{array}

where {\psi }_{1,2} is normalized such that \int {|{\psi }_1|}^{2}dx=1 and \int {|{\psi }_2|}^{2}dx=N_2/N_1. Other parameters are defined as: b_{11}=2a_{11}{N}_1,b_{12}=b_{21}=2m_1{a}_{12}{N}_1/[(1+\gamma )m],b_{22}=2a_{22}{m}_1{N}_1\gamma /m_2,\gamma ={\omega }_{2\perp }/{\omega }_{1\perp },{\lambda }_1={\omega }_{1x}/{\omega }_{1\perp },{\lambda }_2={\omega }_{2x}/{\omega }_{1\perp }, and \kappa =m_1/m_2.

In the absence of external potential and in the case of \kappa =1 and b_{11}=b_{12}=b_{21}=b_{22}, Eq. (14) is perfectly integrable [52]. When integrability is destroyed, closed form can only be given for special cases. In this section, we show that there still exists a specific class of exact vector-soliton solutions for arbitrary interaction strengths, so long as b_{12}^2\ne b_{11}{b}_{22}.

The vector-soliton solution can be found by inserting an appropriate ansatz into Eq. (14). The bright-bright vector soliton solution can be written as:

\begin{array}{}\text{(15)}& \begin{array}{c}{\psi }_{1\mathrm{B}}=\eta \sqrt{C_1}\mathrm{sech}(\eta x-\eta vt){\mathrm{e}}^{\mathrm{i}\left[vx+\left({\eta }^2-v^2\right)t/2\right]}\hfill \\ {\psi }_{2\mathrm{B}}=\eta \sqrt{-C_2}\mathrm{sech}(\eta x-\eta vt){\mathrm{e}}^{\mathrm{i}\left[vx/\kappa +\left(\kappa {\eta }^2-v^2/\kappa \right)t/2\right]}\hfill \end{array}\end{array}

The BD solution is given by

\begin{array}{}\text{(16)}& \begin{array}{c}{\psi }_{1\mathrm{B}}=\eta \sqrt{C_1}\mathrm{sech}(\eta x-\eta vt){\mathrm{e}}^{\mathrm{i}\left[vx+f_{1}t\right]}\hfill \\ {\psi }_{2\mathrm{D}}=[\mathrm{i}v\sqrt{C_2}/\kappa +\eta \sqrt{C_2}\tanh (\eta x-\eta vt)]{\mathrm{e}}^{\mathrm{i}{f}_{2}t}\hfill \end{array}\end{array}

with f_1=\left({\eta }^2-v^2\right)/2-b_{12}{C}_2({\eta }^2+v^2/{\kappa }^2) and f_2=-b_{22}{C}_2({\eta }^2+v^2/{\kappa }^2).

The DD solution is found to be

\begin{array}{}\text{(17)}& \begin{array}{c}{\psi }_{1\mathrm{D}}=[\mathrm{i}v\sqrt{-C_1}+\eta \sqrt{-C_1}\tanh (\eta x-\eta vt)]{\mathrm{e}}^{\mathrm{i}{f}_{1}t}\hfill \\ {\psi }_{2\mathrm{D}}=[\mathrm{i}v\sqrt{C_2}/\kappa +\eta \sqrt{C_2}\tanh (\eta x-\eta vt)]{\mathrm{e}}^{\mathrm{i}{f}_{2}t}\hfill \end{array}\end{array}

with f_1=-{\eta }^2-v^2(-b_{11}{C}_1+b_{12}{C}_2/{\mu }^2) and f_2=-\mu {\eta }^2-v^2(-b_{21}{C}_1+b_{22}{C}_2/{\mu }^2).

In Eqs. (15), (16) and (17), C_1\equiv (b_{22}-\kappa b_{12})/(b_{12}^2-b_{11}{b}_{22}),C_2\equiv (b_{12}-\kappa b_{11})/(b_{12}^2-b_{11}{b}_{22}), the parameter η determines the width of the soliton and can be found by the normalization condition for {\psi }_{1,2}. The parameter v gives the velocity of the soliton. It should be mentioned that the exact vector-soliton solutions are obtained under the assumption that the quantities under a radical sign must be a number greater than zero; thus giving the constrains for the existence of the exact vector-soliton solutions.

To translate our results relevant to the experiments, we take ⁷Li as species 1 and ²³Na as species 2, with fixed and realistic intra-species interaction strength and variable inter-species interaction strength b_{12}. Figure 5 depicts the “phase diagram” for different regimes of vector solitons as b_{12} is varied.

Next, we assume m_1=m_2,{\omega }_1={\omega }_2=\omega, and define variables {\sigma }_1=b_{12}^2-b_{11}{b}_{22},{\sigma }_2=b_{12}-b_{11}, and {\sigma }_3=b_{12}-b_{22} for the facilitation of expression. We first consider the simplest case with the interactions taking exactly the same form. In this case, these coupled equations are nothing but the well-known integrable model proposed by Manakov [53] or the defocusing–defocusing NLS equations. Then we can directly use the method developed in subsection 2.2, and obtain a similar integrability condition for the existence of the exact symmetric vector solitons [48].

We now turn to Eq. (14), with unequal interaction parameters b_{11}(t)\ne b_{22}(t)\ne b_{12}(t). We take advantage of the extend mapping deformation method, which does not depend on the integrability, to obtain the exact vector-soliton solutions for the general Eq. (14).

Expressing the order parameters in terms of their modulus and phases, i.e., {\psi }_i(x,t)=\sqrt{n_i}\exp (\mathrm{i}{\Delta }_i), and then separating real and imaginary parts, we obtain a set of coupled nonlinear equations for n_i and {\Delta }_i

\begin{array}{}\text{(18)}& \begin{array}{c}{\sqrt{n_i}}_{xx}-\sqrt{n_i}{\Delta }_{it}-\sqrt{n_i}{\Delta }_{ix}^2-b_{12}{n}_{(3-i)}\sqrt{n_i}-b_{ii}{(\sqrt{n_i})}^3-{\lambda }^2{x}^2\sqrt{n_i}=0\hfill \\ {\sqrt{n_i}}_t+\sqrt{n_i}{\Delta }_{ixx}+2{\sqrt{n_i}}_x{\Delta }_{ix}=0\hfill \end{array}\end{array}

where i=1,2. By setting \sqrt{n_i}=A_0(t)+A_i(t)\varphi (\xi ),\xi =p_1(t)x+p_2(t), and {\Delta }_i=k_i(t)+{\Gamma }_i(t)x+{\Lambda }_i(t)x^2 and using the auxiliary equation {\varphi }^{{}_{}{}^{\prime }2}=c_0+c_2{\varphi }^2+c_4{\varphi }^4, we derive a set of over-determined partial differential equations with respect to p_1(t),p_2(t),k_i(t),{\Gamma }_i(t),{\Lambda }_i(t),A_0(t),A_1(t), and A_2(t). Finally, solving these equations, we can obtain the following constraints for the existence of the exact vector-soliton solutions:

\begin{array}{}\text{(19)}& \begin{array}{c}{b}_{11}{A}_1^2(t)+b_{12}{A}_2^2(t)=b_{22}{A}_2^2(t)+b_{12}{A}_1^2(t)\hfill \\ \frac{-{\partial }^2{p}_1(t)}{\partial t^2}{p}_1(t)+2{\left[\frac{\partial p_1(t)}{\partial t}\right]}^2+4p_1^2(t){\lambda }^2=0\hfill \\ p_1^2(t)-b_{11}{p}_1(t)-b_{12}\frac{{\sigma }_2}{{\sigma }_3}{p}_1(t)=0\hfill \end{array}\end{array}

Then, we get a relation between the the scattering length and trapping potential, which reads

\begin{array}{}\text{(20)}& \frac{b_{12}-b_{22}}{b_{12}^2-b_{11}{b}_{22}}=\frac{a_1\sin (2\lambda t)-a_2\cos (2\lambda t)}{2\lambda }\end{array}

In the case of all the interaction parameters , Eq. (14) has bright–bright vector solitons

\begin{array}{}\text{(21)}& \begin{array}{c}{\psi }_{1\mathrm{B}}=C_1\sqrt{-a(t)}\mathrm{sech}(\xi (x,t))\exp \left(\mathrm{i}{v}_{1\mathrm{B}}(x,t)\right)\hfill \\ {\psi }_{2\mathrm{B}}=C_2\sqrt{-\sigma a(t)}\mathrm{sech}(\xi (x,t))\exp \left(\mathrm{i}{v}_{2\mathrm{B}}(x,t)\right)\hfill \end{array}\end{array}

where

When a(t)>0 and \sigma (t)>0, Eq. (14) has dark–dark vector solitons

\begin{array}{}\text{(23)}& \begin{array}{c}{\psi }_{1\mathrm{D}}=D_1\sqrt{a(t)}\tanh (\xi (x,t))\exp \left(\mathrm{i}{v}_{1\mathrm{D}}(x,t)\right)\hfill \\ {\psi }_{2\mathrm{D}}=D_2\sqrt{\sigma a(t)}\tanh (\xi (x,t))\exp \left(\mathrm{i}{v}_{2\mathrm{D}}(x,t)\right)\hfill \end{array}\end{array}

Thus, we conclude the additional conditions for the existence of the exact vector-soliton solutions

with {\Delta }_1\equiv -{\sigma }_3/{\sigma }_1 and {\Delta }_2\equiv {\sigma }_2/{\sigma }_1.

Shown in Table 1 are the regions where the dark-dark vector solitons can exist when the intercomponent interaction parameter b_{ij} is varied and intracomponent ones are constants. Case 1 has been extensively studied in Refs. [44, 54]. It is interesting to observe that in cases 2 and 3, a dark soliton can be formed in a condensate with self-attractive interaction. This can be understood by the fact that in the two-component BECs, the repulsive interaction coming from the second component or the intercomponent interaction induces an effective repulsive interaction in the self-attractive one, which leads to the formation of dark solitons in the component with self-attractive interaction. In addition, we would like to point out that the formation of this type of vector solitons has a threshold value in b_{12}, which indicates that the intercomponent repulsive interaction must be strong to induce an effective repulsive interaction in the self-attractive ones [55].

To study the dynamic stability of the vector-soliton solutions [Eqs. (15), (16) and (17)], we calculate the collective excitation frequencies by solving the corresponding Bogoliubov-de Gennes equation [56]. The results show that the excitation frequencies of the BB and BD vector-solitons are all real, while those of the DD vector-soliton contain complex values. This stability analysis is confirmed by direct numerical simulations of Eq. (14). Computational examples are shown in Fig. 6. To make the simulations more realistic, we add a harmonic trapping potential along the longitudinal direction. As to the initial condition, we use the exact solution for the bright component, while multiplying Eq. (16) by a Thomas–Fermi profile to simulate the dark component. The upper panel of Fig. 6 shows the evolution of a BD vector-soliton. The presence of the weak longitudinal trap causes the vector-soliton to oscillate while maintaining the overall shape. By contrast, the lower panel of Fig. 6 shows that the DD soliton is dynamically unstable — an initial DD vector-soliton is quickly destroyed and the two species tend to phase separate.

Figure 7 illustrates the conversion between different types of vector-solitons via Feshbach resonance. The initial condition is a zero-velocity BD vector-soliton. The value of b_{12} is then decreased and the system evolves to two BB vector-soliton pairs. Then, we restore the initial value of b_{12} before the two BB vector-solitons merge into one. As seen in Fig. 7, such a controlled tuning of b_{12} can also recover the initial BD vector-soliton. This represents a remarkable example of the control over the BECs dynamics, afforded by the tunability of inter-species interaction strength [57].

2.5 Dynamics of two-component BECs with time and space-modulated nonlinearities

Now, we turn to the two-component BECs with time and space-modulated nonlinearities which is governed by the coupled one-dimensional GP equations [57-59]:

\begin{array}{}\text{(25)}& \begin{array}{c}\mathrm{i}\frac{\partial {\psi }_1}{\partial t}=(-\frac{1}{2}\frac{{\partial }^2}{\partial x^2}+\frac{{\gamma }^2}{2}{x}^2+b_{11}{|{\psi }_1|}^2+b_{12}{|{\psi }_2|}^2){\psi }_1\hfill \\ \mathrm{i}\frac{\partial {\psi }_2}{\partial t}=(-\frac{1}{2}\frac{{\partial }^2}{\partial x^2}+\frac{{\gamma }^2}{2}{x}^2+b_{12}{|{\psi }_1|}^2+b_{22}{|{\psi }_2|}^2){\psi }_2\hfill \end{array}\end{array}

where b_{11}=2a_{11},b_{22}=2a_{22},b_{12}=2a_{12} with the scattering lengths a_{11},a_{22} and a_{12} being spatiotemporally inhomogeneous, which are experimentally feasible due to the flexible and precise control of the scattering lengths achievable in quasi-1D BECs with magnetically tuning the Feshbach resonances.

We consider the exact spatially localized solutions of Eq. (25) for which {lim}_{|x|\to \mathrm{\infty }}{\psi }_{1,2}(x,t)=0, taking the similarity transformation

\begin{array}{}\text{(26)}& {\psi }_i(x,t)={\beta }_i(x,t){\mathrm{e}}^{\mathrm{i}{\alpha }_i(x,t)}U[X(x,t)],i=1,2\end{array}

with U(X),V(X) satisfying

\begin{array}{}\text{(27)}& \begin{array}{c}{U}_{XX}+g_{11}{U}^3+g_{12}U{V}^2=0\hfill \\ V_{XX}+g_{22}{V}^3+g_{12}V{U}^2=0\hfill \end{array}\end{array}

As aforementioned, we define {\sigma }_1=g_{12}-g_{11},{\sigma }_2=g_{12}-g_{22} and {\sigma }_{12}=g_{12}^2-g_{11}{g}_{22}. Substituting Eq. (26) into Eq. (25) and asking U(X),V(X) to satisfy Eq. (27) we conclude that when

\begin{array}{}\text{(28)}& b_{11}=g_{11}\theta (t,x),b_{12}=g_{12}\theta (t,x),b_{22}=g_{22}\theta (t,x)\end{array}

with \theta (t,x)=-{\lambda }^2{\mathrm{e}}^{-3{\lambda }^2{x}^2-6\lambda \delta x-2{\delta }^2}/[2{\zeta }_3^2\left(t\right)], we can obtain the following conditions (for i=1,2)

\begin{array}{}\text{(29)}& \begin{array}{c}{\alpha }_i={\zeta }_i\left(t\right)-{\lambda }_t{x}^2/2\lambda -{\delta }_{t}x/\lambda \hfill \\ {\beta }_1={\beta }_2={\zeta }_3\left(t\right){\mathrm{e}}^{\lambda x(\lambda x+2\delta )/2}\hfill \\ X=\sqrt{\pi }\mathrm{erf}\left(\lambda x+\delta \right)/2\hfill \end{array}\end{array}

with \mathrm{erf}(s)=\frac{2}{\sqrt{\pi }}{\int }_0^s{\mathrm{e}}^{-{\tau }^2}d\tau ,{\zeta }_i\left(t\right) and {\zeta }_3\left(t\right) satisfy

\begin{array}{}\text{(29)}& \begin{array}{c}{\zeta }_i\left(t\right)={\displaystyle \int ({\lambda }^4-\delta _t^{}{}_{}{}^2+{\lambda }^4{\delta }^2)/(2{\lambda }^2)dt}+C_i\hfill \\ {\zeta }_3\left(t\right)=\sqrt{\lambda }{\mathrm{e}}^{\frac{1}{2}{\delta }^2}\hfill \end{array}\end{array}

with C₁ and C₂ arbitrary constants and

\begin{array}{}\text{(30)}& \begin{array}{c}{\lambda }^6-{\gamma }^2{\lambda }^2-2\lambda _t^{}{}_{}{}^2+{\lambda }_{tt}\lambda =0\hfill \\ 2{\delta }_{tt}\lambda +2{\lambda }^5\delta -4{\delta }_t{\lambda }_t=0\hfill \end{array}\end{array}

When setting {\sigma }_1/{\sigma }_{12}>0 and {\sigma }_2/{\sigma }_{12}>0, Eq. (27) has two families of exact solutions as:

\begin{array}{}\text{(31)}& \begin{array}{c}{U}^{(1)}(X)=\sqrt{{\sigma }_2/{\sigma }_{12}}{\nu }_1\mathrm{cn}({\nu }_{1}X,\sqrt{2}/2)\hfill \\ V^{(1)}(X)=\sqrt{{\sigma }_1/{\sigma }_{12}}{\nu }_1\mathrm{cn}({\nu }_{1}X,\sqrt{2}/2)\hfill \end{array}\end{array}

and

\begin{array}{}\text{(32)}& \begin{array}{c}{U}^{(2)}(X)=\sqrt{2}/2\sqrt{{\sigma }_2/{\sigma }_{12}}{\nu }_2\mathrm{sd}({\nu }_{2}X,\sqrt{2}/2)\hfill \\ V^{(2)}(X)=\sqrt{2}/2\sqrt{{\sigma }_1/{\sigma }_{12}}{\nu }_2\mathrm{sd}({\nu }_{2}X,\sqrt{2}/2)\hfill \end{array}\end{array}

where {\nu }_1,{\nu }_2 are arbitrary constants. When imposing the bounded condition {lim}_{|x|\to \mathrm{\infty }}{\psi }_{1,2}(x,t)=0, we have {\nu }_1=2(2n+1)K\left(\frac{\sqrt{2}}{2}\right)/\sqrt{\pi } for Eq. (31) and {\nu }_2=4mK\left(\frac{\sqrt{2}}{2}\right)/\sqrt{\pi } for Eq. (32), where n and m are integer numbers and K\left(\frac{\sqrt{2}}{2}\right)={\int }_0^{\pi /2}{[1-{\left(\frac{\sqrt{2}}{2}\right)}^2{\sin }^2\tau ]}^{-1/2}d\tau is elliptic integral of the first kind.

From Eq. (30), we have

\begin{array}{}\text{(33)}& \delta \left(t\right)=c_1{\mathrm{e}}^{\mathrm{i}\int {\lambda }^{2}dt}+c_2{\mathrm{e}}^{-\mathrm{i}\int {\lambda }^{2}dt}\end{array}

with c₁ and c₂ arbitrary constants.

If we set \lambda =1/\xi , the first one of Eq. (30) becomes

\begin{array}{}\text{(34)}& {\xi }_{tt}+{\gamma }^2\xi =1/{\xi }^3\end{array}

To obtain the solutions of Eq. (34) we let γ satisfy

\begin{array}{}\text{(35)}& {\gamma }^2={\gamma }_0^2+ϵ\cos ({\gamma }_{1}t)\end{array}

Therefore, the general solution to Eq. (30) is

\begin{array}{}\text{(36)}& \lambda ={(A{\xi }_1^2+B{\xi }_2^2+2C{\xi }_1{\xi }_2)}^{-1/2}\end{array}

where A,B,C are constants satisfying AB-C^2=1/W^2, and the Wronskian W={\xi }_1{\xi }_{2t}-{\xi }_2{\xi }_{1t} with ({\xi }_1,{\xi }_2) being two linearly independent solutions of homogeneous ordinary differential equation

\begin{array}{}\text{(37)}& {\xi }_{tt}+{\gamma }^2\xi =0\end{array}

Combining Eqs. (26) and (29) with (31) and (32), we arrive at two families of explicitly exact solutions for the coupled 1D GP Eqs. (25) as:

\begin{array}{}\text{(38)}& \begin{array}{c}{\psi }_1^{(j)}(x,t)=\sqrt{\lambda }{\mathrm{e}}^{\frac{1}{2}{\delta }^2}{\mathrm{e}}^{\frac{\lambda x(\lambda x+2\delta )}{2}}{\mathrm{e}}^{\mathrm{i}{\alpha }_1(x,t)}{U}^{(j)}(X)\hfill \\ {\psi }_2^{(j)}(x,t)=\sqrt{\lambda }{\mathrm{e}}^{\frac{1}{2}{\delta }^2}{\mathrm{e}}^{\frac{\lambda x(\lambda x+2\delta )}{2}}{\mathrm{e}}^{\mathrm{i}{\alpha }_2(x,t)}{V}^{(j)}(X)\hfill \end{array}\end{array}

where U^{(j)}(X) and V^{(j)}(X) are given by Eqs. (31) and (32), index j=1,2,{\alpha }_1,{\alpha }_2 satisfy Eq. (29), and X,\delta ,\lambda satisfy Eqs. (29), (33) and (36), respectively.

In order to investigate the dynamics of the explicitly exact solutions (38), we take special parameters {\gamma }_0,{\gamma }_1 and ϵ. When the ratio of the confining frequencies \gamma is time-independent, that is, parameters {\gamma }_1=ϵ=0, solving Eq. (37) we have \lambda =1/\xi with \xi ={[A-(A-B){\cos }^2({\gamma }_{0}t)+\sqrt{AB\gamma _0^{}{}_{}{}^2-1}\sin (2{\gamma }_{0}t)/{\gamma }_0]}^{1/2}. Here \xi is the width of the explicitly exact solutions (38) and \sqrt{\lambda } is its amplitude. Specially, we further suppose c_1=c_2=0, i.e., \delta =0. When the cigar-shaped trap with axial frequency {\omega }_x=70\pi Hz and radial frequency {\omega }_{\perp }=800\pi Hz is considered, the ratio of the confining frequencies \gamma ={\omega }_x/{\omega }_{\perp }=7/80. Therefore, from Eq. (35) the parameter {\gamma }_0=\gamma =7/80. In order to determine parameters A and B in Eq. (36), we further consider the initial condition of Eq. (34) as \xi (0)=\sqrt{3} and {\xi }_t(0)=1/80. In the following, within the safe region we will always take the parameters g_{ij} in Eq. (28) to be g_{11}=-1,g_{22}=-3,g_{12}=2.

In Fig. 8, we show the evolution of condensate density of the wave functions {\psi }_1^{(1)} and {\psi }_1^{(2)}. Figure 8(a) and (b) demonstrate the density profiles of the wave function {\psi }_1^{(1)} for n=0,1, respectively, and Fig. 8(c) and (d) demonstrate the density profiles of the wave function {\psi }_1^{(2)} for m=1,2, respectively. Figure 8(e) demonstrates the width and amplitude of the wave functions. It is observed that the localized nonlinear matter waves are space localized and time periodic, which are usually called breathing solitons. Here, n and m are the order of the breathing solitons. It is also observed that the amplitude and width of the localized nonlinear waves vary periodically about time.

When \gamma is a function of time, we can select proper parameters {\gamma }_0,{\gamma }_1 and ϵ to formulate a type of quasi-breathing solitons. Specially, we choose parameters {\gamma }_0=7/80, {\gamma }_1=\sqrt{2}, and ϵ=1/2. In this case, the two solutions of Eq. (37) are Mathieu functions {\xi }_1=\mathrm{MathieuC}(49/3200,-1/2,\sqrt{2}t/2),{\xi }_2=\mathrm{MathieuS}(49/3200,-1/2,\sqrt{2}t/2) which are in the stability region of Eq. (37). We will consider the initial condition of Eq. (34) to be \xi (0)=\sqrt{3} and {\xi }_t(0)=\sqrt{6}/6. In Figure 9, we show the development of density profiles for the wave functions {\psi }_1^{(1)} and {\psi }_2^{(2)} in (38). Figure 9(a) and (b) demonstrate the density profiles of the wave functions {\psi }_1^{(1)} for order n=0, and Fig. 9(c) and (d) demonstrate the density profiles of the wave functions {\psi }_2^{(2)} for order m=1. It is observed that the localized nonlinear matter waves are space localized and time quasi-periodic, which are quasi-breathing solitons with orders n or m.

Next, we let the parameters \delta =0,{\gamma }_0=7/80, {\gamma }_1=4 and ϵ=1/2. In this case, we are in the instability region of Eq. (34). However, the solutions of Eq. (37) are also Mathieu functions {\xi }_1=\mathrm{MathieuC}(49/25600,-1/16,2t),{\xi }_2=\mathrm{MathieuS}(49/25600,-1/16,2t). After some algebra, we also get the analytical expression of function \lambda . In order to determine the unknown parameters A and B in λ, we choose the initial condition of Eq. (34) as \xi (0)=1 and {\xi }_t(0)=1. Figure 10(a) and (b), and (c) and (d) describe the evolution of the condensate density profiles for the wave functions {\psi }_1^{(1)} with order n=0 and {\psi }_2^{(2)} with order m=1, respectively. It is shown that the density wave packets are space localized and time resonant, which have the behaviors of resonant solitons. At the beginning, resonant solitons are localized nonlinear matter waves with low amplitude and large width, after some time, their amplitudes become high but widths become small; as time goes on their amplitudes decrease and widths increase gradually. The nonlinear matter waves demonstrate transitory resonant soliton behaviors.

Then we consider the case of \delta \ne 0. In this case, the nonlinearities in Eq. (25) become more complicated, the amplitude of the nonlinear matter wave becomes \sqrt{\lambda }{\mathrm{e}}^{\frac{1}{2}{\delta }^2}, and the center of the solitons can move following the time because of variable X=\frac{\sqrt{\pi }}{2}\mathrm{erf}\left(\lambda x+\delta \right). So for \delta \ne 0, our localized nonlinear matter waves are localized moving solitons. To do so, we still choose {\gamma }_0=7/80 and {\gamma }_1=ϵ=0. Figure 11(a) and (b) describe the evolutions of the density profiles for the wave function {\psi }_1^{(1)} with orders n=0,1, respectively; Fig 11(c) and (d) describe the evolutions of the density profiles for the wave function {\psi }_1^{(2)} with orders m=1,2, respectively. Figure 11(e) demonstrates the shapes of the function \sqrt{\lambda }, the width \xi =1/\lambda and the amplitude \sqrt{\lambda }{\mathrm{e}}^{\frac{1}{2}{\delta }^2} of the moving breathing solitons. It is found that the localized nonlinear matter waves are space localized and time periodically moving and the amplitudes of the moving breathing solitons are higher than the breathing solitons in Fig. 8.

Finally, we study the dynamical stability of the exact localized nonlinear wave solutions (38) for other parameters by adding certain initial stochastic noises with various intensities such as Gaussian distributed noise and uniform distribution noise. Our numerical calculations show that only for order n=0 (ground state) is the breathing soliton solution (38) with \delta =0 dynamically stable, while the moving soliton solution (38) with \delta \ne 0 is dynamically unstable for all n [58].

3 Dynamics of BECs in two dimension

3.1 Dynamics of ring dark soliton in BECs with tunable interaction

We consider a BECs trapped in an external harmonic potential V(\mathit{r},z)=m({\omega }_r^2{r}^2+{\omega }_z^2{z}^2)/2, where r^2=x^2+y^2. When \Omega \equiv {\omega }_r/{\omega }_z\ll 1, the system is governed by the quasi-2D GP equation

\begin{array}{}\text{(39)}& \mathrm{i}{\displaystyle \frac{\partial \psi }{\partial t}}=-\frac{1}{2}{\nabla }^2\psi +g(t){\mid \psi \mid }^2\psi +\frac{1}{2}{\Omega }^2{r}^2\psi \end{array}

where the length, time, and wave function are scaled as x=a_h\tilde{x}, t=\frac{\tilde{t}}{{\omega }_z}, \psi =\frac{\psi }{a_h\sqrt{4\pi a_0\eta }}. {\nabla }^2={\partial }^2/\partial x^2+{\partial }^2/\partial y^2={\partial }^2/\partial r^2+1/r\times \partial /\partial r+{\partial }^2/\partial {\theta }^2, \Omega ={\omega }_r/{\omega }_z, g(t)=a(t)/a_0.

We first consider the analytic solution of Eq. (39) with circular symmetry, \psi (r,t). When there is no trap \Omega =0 and the interaction is time-independent, g(t)=C, the system is governed by the standard NLS equation:

\begin{array}{}\text{(40)}& \mathrm{i}\frac{\mathrm{d}Q(R,T)}{\mathrm{d}T}+\frac{1}{2}[\frac{{\partial }^{2}Q(R,T)}{\partial R^2}+\frac{1}{R}\frac{\partial Q(R,T)}{\partial R}]-C{\mid Q(R,T)\mid }^{2}Q(R,T)=0\end{array}

Under small-amplitude approximation, Eq. (40) can be transformed to the famous cylindrical KdV (cKdV) equation, whose exact solutions are well studied.

For the more general case, especially when there is one time-dependent trap V(t), the transformation method has to be used. By using a transformation: \psi (r,t)=Q(R(r,t),T(t)){\mathrm{e}}^{\mathrm{i}a(r,t)+c(t)}, the general form of Eq. (39) can be converted to Eq. (40). In the transformation, R(r,t), T(t), a(r,t), and c(t) are assumed to be real functions and read R(r,t)=\alpha (t)r,T(t)=\int {\alpha }^2(t^{\prime })d{t}^{\prime }+C_0, c(t)=\frac{1}{2}\ln \frac{{\alpha }^2(t)}{C}, a(r,t)=-\frac{1}{2\alpha (t)}\frac{\mathrm{d}\alpha (t)}{\mathrm{d}t}{r}^2. Meanwhile, \alpha (t) and \Omega satisfy the following condition:

\begin{array}{}\text{(41)}& {\displaystyle \frac{1}{\alpha (t)}}{\displaystyle \frac{{\mathrm{d}}^2\alpha (t)}{\mathrm{d}{t}^2}}-\frac{2}{\alpha {(t)}^2}{\left(\frac{\mathrm{d}\alpha (t)}{\mathrm{d}t}\right)}^2-{\Omega }^2=0\end{array}

If we substitute y(t)=\mathrm{d}\alpha /(\alpha \mathrm{d}t) into Eq. (41), the equation becomes \frac{\mathrm{d}y(t)}{\mathrm{d}t}=y^2+{\Omega }^2. This equation is nothing but the standard Riccati equation. It can be solved not only for \Omega =constant, but also for various types of {\Omega }^2(t), such as a+b\sin (\lambda t), a+b\cosh (t), a{\mathrm{e}}^{\lambda t} (a, b, λ are arbitrary constants) and so on. Thus, we can study the dynamics of a system with time-dependent external trap analytically.

When the RDSs get deeper, the small-amplitude approximation is invalid and we have to appeal to the numerical simulation. In the following sections, the dynamics and stability of RDS is studied numerically. We first consider a simple case with g(t)=1, \Omega =0.028, R_0=28.9. These parameters corresponding the parameters in real experiment: {}_{}{}^{87}\mathrm{Rb} (a_0=5.7 nm) condensate of radius 30 μm, containing 20000 atoms in a disk-shaped trap with {\omega }_r=2\pi \times 18 Hz and {\omega }_z=2\pi \times 628 Hz. Also, the RDS has the radius R₀ = 12.38 μm.

Our numerical results show that the shallow RDS is stable against the small distortion of the ringshape. Furthermore, we find that there exists a maximum eccentricity e_{\max } for a given initial depth. When , the RDSs are stable and oscillate reserving their shape with the same period as the unperturbed RDSs until decaying out. With the increasing of the depth, the e_{\max } becomes smaller. On the contrary, when e_c exceeds e_{\max }, snaking sets in and the dark soliton breaks into two vortex pairs, presenting a striking contrast to the multiples of four pairs reported previously. Thus, we take the following equation as the initial configuration:

\begin{array}{}\text{(42)}& \psi (x,y,0)=(1-{\Omega }^2{r}^2/4)\cdot \left[\cos \varphi (0)\tanh Z(r_1)+\mathrm{i}\sin \varphi (0)\right]\end{array}

where Z(r_1)=(r_1-R_0)\cos \varphi (0), r_1=\sqrt{(1-e_c^2)x^2+y^2}, e_c is the eccentricity of the ring and \cos \varphi (0) is proportional to the depth of the input soliton. When e_c\ne 0, R₀ represents the length of semiminor axis of the elliptical configuration.

To illustrate our results, we consider one typical case with \cos \varphi (0)=0.6 and e_c=0.4. As shown in Figs. 12 and 13, the ring dark soliton initially shrinks and when reaching the minimum radius in the short-axis direction, it starts snaking and forming two dark lumplike solitons in the horizontal direction; they move in the opposite direction and then break into two vortex pairs. Simultaneously the vortices and antivortices move along the ring. The result of motion is their collision in pairs in the vertical direction, followed by vortices and antivortices’s mergence and totally forming a pair of lumplike solitons. Then the lumplike solitons move towards to each other, (to the center of condensate), in the vertical direction. When they reach the minimum radius, a ring dark soliton forms. After a short time the system returns to the state of two dark lumplike solitons. The lumplike solitons leave each other and in the path, each of the lumplike solitons breaks into a vortex pair again. However, the configuration is different to the initial one: the vortex (antivortex) is substituted by antivortex (vortex), just like the vortex and the antivortex passing each other directly following the motion before the mergence. Then the vortices and the antivortices keep on moving, then collide and merge in the horizontal direction. The dynamical process repeats itself.

Furthermore, we found that the Feshbach resonance can dramatically affect the dynamics of dark solitons. We choose e_c=0, and g(t)\sim a{\mathrm{e}}^{\pm \omega t} or g(t)\sim a\pm \sin (\omega t), where a and ω are arbitrary constant parameters. Our results show that for the deep ring dark solitons, the Feshbach resonance management remarkably changes the evolution and the instability of solitons. Shown in Fig. 14(a), the black ring dark solitons break into 20 vortex pairs; while in Fig. 14(b), there are only 4 pairs. What is more interesting and more important are the cases in which the lifetime of the RDS can be extended, such as g(t)=1-\sin (\omega t), as shown in Table 2, where the lifetime of RDSs corresponding to different modulation frequency of the FRM (ω) is given [60].

3.2 Quantized quasi-two-dimensional Bose–Einstein condensates

We also consider a BECs confined in a harmonic trap, but with spatially modulated nonlinearity. In this case, the GP equation reads as:

\begin{array}{}\text{(43)}& \mathrm{i}{\psi }_t=-\frac{1}{2}({\psi }_{xx}+{\psi }_{yy})+\frac{1}{2}{\omega }^2(x^2+y^2)\psi +g(x,y){|\psi |}^2\psi \end{array}

where g(x,y)=4\pi a_s(x,y) can be spatially inhomogeneous by magnetically tuning the Feshbach resonances.

Now, we consider the spatially localized stationary solution \psi (x,y,t)=\varphi (x,y){\mathrm{e}}^{-\mathrm{i}\mu t} of Eq. (43) with \varphi (x,y) being a real function for {lim}_{|x|,|y|\to \mathrm{\infty }}\varphi (x,y)=0. Solving the stationary equation by similarity transformation [61], we obtain a family of exact localized nonlinear wave solutions for Eq. (43) as:

\begin{array}{}\text{(44)}& {\psi }_n={\displaystyle \frac{(n+1)K(k)\eta }{\sqrt{\nu }}}\mathrm{cn}(\theta ,k){\mathrm{e}}^{-\mathrm{i}\mu t},n=0,2,4,\mathrm{\cdots }\end{array}

\begin{array}{}\text{(45)}& {\psi }_n={\displaystyle \frac{(n+1)K(k)\eta }{\sqrt{2\nu }}}\mathrm{sd}(\theta ,k){\mathrm{e}}^{-\mathrm{i}\mu t},n=1,3,5,\mathrm{\cdots }\end{array}

where k=\sqrt{2}/2 is the modulus of elliptic function, ν is a positive real constant, K(k)={\int }_0^{\frac{\pi }{2}}{(1-k^2{\sin }^2\varsigma )}^{-1/2}d\varsigma is elliptic integral of the first kind, \theta ,\eta and g are determined by

\begin{array}{}\text{(46)}& \begin{array}{c}\theta =(n+1)K(k)\mathrm{erf}[\sqrt{2\omega }\left(x+y\right)/2]\hfill \\ \eta =e^{\omega xy}\mathrm{KummerU}[-\mu /(2\omega ),1/2,\omega {\left(x-y\right)}^2/2]\hfill \\ g(x,y)=-2\omega \nu /(\pi {\eta }^2){\mathrm{e}}^{-\omega {\left(x+y\right)}^2}\hfill \end{array}\end{array}

here \mathrm{erf}(x)=\frac{2}{\sqrt{\pi }}{\int }_0^x{\mathrm{e}}^{-{\tau }^2}d\tau is error function, and KummerU(a,c,s) [62] is Kummer function of the second kind which is a solution of ordinary differential equation s{\Lambda }^{\prime \prime }(s)+(c-s){\Lambda }^{\prime }(s)-a\Lambda (s)=0. It is easy to see that when |x|,|y|\to \mathrm{\infty } we have {\psi }_n\to 0 for solutions {\psi }_n in Eqs. (44) and (45) with Eq. (46), thus they are localized bound state solutions.

From Eq. (46), we observe that the number of zero points of function η is equal to that of function \mathrm{KummerU}[-\mu /(2\omega ),1/2,\omega {\left(x-y\right)}^2/2]. We assume the number of zero points in η along line y=-x is l. In the following, we will see that integer n is associated with the energy levels of the atoms and integers n,l determine the topological properties of atom packets, so n and l are named the principal quantum number and secondary quantum number in quantum mechanics. In addition, the three free parameters ω, μ and ν are positive, so the dimensionless interaction function g(x,y) is negative, which indicates an attractive interaction between atoms.

Figure 15 shows the density distribution of the quasi-2D BECs with spatially modulated nonlinearities in harmonic potential for l=0. It is easy to see that the matter wave functions [Eq. (44)] satisfy {\psi }_n(-x,-y)={\psi }_n(x,y), so they are even parity and are invariant under space inversion. Figure 15(a)–(c) demonstrate the density profiles of the even parity wave functions [Eq. (44) with Eq. (46)] for n=0,2,4, which correspond to a low energy state and two highly excited states. The matter wave functions [Eq. (45)] satisfy {\psi }_n(-x,-y)=-{\psi }_n(x,y), which denotes that they are odd parity. Figure 15(d)–(f) demonstrate the density profiles of the odd parity wave functions [Eq. (45) with Eq. (46)] for n=1,3,5, which correspond to three highly excited states. It is observed that when the secondary quantum number l=0, the number of nodes along line y=x for each quantum state is equal to the corresponding principal quantum number n. Also, the number of density packets increases one by one along line y=x when the n increases. This is similar to the quantum properties in the linear harmonic oscillator.

Next, when the principal quantum number n is fixed, we tune the secondary quantum number l to observe the novel quantum phenomenon in quasi-2D BECs. In Fig. 16 we demonstrate the density distributions of quasi-2D BECs for different secondary quantum number. Figure 16(a)–(d) show the density profiles of the even parity wave function [Eq. (44) with Eq. (46)] for n=0, and l=0,1,2 and 3, respectively. It is seen that the number of nodes for the density packets along line y=-x is equal to the corresponding secondary quantum number l which describes the topological patterns of the atom packets, and the number of density packets increases one by one when l increases. Figure 16(e)–(h) show the density profiles of the odd parity wave function [Eq. (45) with Eq. (46)] for n=1 and l=0,1,2,3. We see that the number of density packets increases pair by pair when l increases. The number of density packets for each quantum state is equal to (n+1)\times (l+1), and all the density packets are symmetrical with respect to lines y=\pm x, as shown in Figs. 15 and 16 [63].

3.3 Quantized vortex solitons in quasi-two-dimensional Bose–Einstein condensates

The scaled form of the NLSE is

\begin{array}{}\text{(46)}& \mathrm{i}{\psi }_t=-{\nabla }^2\psi +g(r){|\psi |}^2\psi +V(r)\psi \end{array}

where {\nabla }^2 is the 2D Laplacian. Assuming the stationary wave function as \psi (r,\theta ,t)=\varphi (r)\exp (\mathrm{i}S\theta -\mathrm{i}\mu t), where θ is the azimuthal angle, S an integer vorticity, and μ the chemical potential, then we can obtain the equation for the real stationary wave function \varphi (r):

\begin{array}{}\text{(46)}& \mu \varphi =-{\varphi }^{\prime \prime }-r^{-1}{\varphi }^{\prime }+g(r){\varphi }^3+\left[S^2{r}^{-2}+V(r)\right]\varphi \end{array}

For S\ne 0, \varphi (r) should vary as r^{|S|} at r\to 0, which is replaced by {\varphi }^{\prime }(r=0)=0 for S=0. The localization requires \varphi (r=\mathrm{\infty })=0.

After defining \varphi (r)\equiv \rho (r)U(R(r)),g(r)\equiv g_0{r}^{-2}\cdot {\rho }^{-6}(r), with R(r)\equiv {\int }_0^r{s}^{-1}{\rho }^{-2}(s)ds, we can obtain

\begin{array}{}\text{(47)}& {\rho }^{\prime \prime }+r^{-1}{\rho }^{\prime }+\left(\mu -V(r)-S^2{r}^{-2}\right)\rho =E{r}^{-2}{\rho }^{-3}\end{array}

\begin{array}{}\text{(48)}& -{\mathrm{d}}^{2}U/\mathrm{d}{R}^2+g_0{U}^3=EU\end{array}

where E and g₀ are constants. Then, we can construct exact solutions to the underlying GP equation if a solution to Eq. (47) is known, and the nonlinearity-modulation profile admitting exact solutions is determined by \rho (r). From the view point of physics, ρ cannot change its sign; further, it must behave as r^{-a} with a⩾1/3 at r\to 0, and diverge (\rho \to \mathrm{\infty }) at r\to \mathrm{\infty }, so that the nonlinearity strength is bounded and the integration in R(r) converges.

When E=0, one can construct the exact VS solutions for the attractive nonlinearity (), and Eq. (47) is solvable. In the case of harmonic potential V=k{r}^2, ρ can be found in terms of the Whittaker’s M and W functions [64]: \rho (r)=r^{-1}[c_1\mathrm{M}(\mu /4\sqrt{k},|S|/2,\sqrt{k}{r}^2)+c_2\mathrm{W}(\mu /4\sqrt{k},|S|/2,\sqrt{k}{r}^2)], where the restrictions on ρ require and c_1{c}_2>0. On the contrary, without the potential,(k=0), ρ degenerates to \rho (r)=c_3{I}_S(\sqrt{-\mu }r)+c_4{K}_S(\sqrt{-\mu }r), with , I_S and K_S being the modified Bessel and Hankel functions, and the constants satisfying c_3{c}_4>0. In both cases, one has \rho (r)\sim r^{-|S|} at r\to 0, hence g(r)\sim r^{6|S|-2} and R(r)\sim r^{2|S|} at r\to 0, and \rho (r)\to \mathrm{\infty } as r\to \mathrm{\infty }. Thus, the respective nonlinearity is localized and bounded, and R(r) is bounded too. To meet boundary conditions \varphi (0)=\varphi (\mathrm{\infty })=0, an exact solution to Eq. (48) is chosen as

\begin{array}{}\text{(49)}& U(R)=\left(n\eta /\sqrt{-g_0}\right)\mathrm{cn}(n\eta R-K(1/\sqrt{2}),1/\sqrt{2})\end{array}

where n=2,4,6,\mathrm{\cdots } with n/2 being the radial quantum number, \eta \equiv K(\sqrt{2}/2)/R(r=\mathrm{\infty }), and K(1/\sqrt{2}) is the complete elliptic integral of the first kind.

From Eq. (49), it is easy to see that U(R)\sim R at R\to 0, which implies that the amplitude of the exact VS is \rho U\sim r^{|S|} at r\to 0. Thus, for given μ, S, and nonlinearity strength g₀ (and k, in the presence of the trap), one can construct an infinite number of exact VSs with n/2 bright rings surrounding the vortex core, as shown in Fig. 17.

In order to investigate the linear stability of the obtained solutions, we consider perturbed solutions in the form of \psi (r,\theta ,t)=\exp (\mathrm{i}S\theta -\mathrm{i}\mu t)[\varphi (r)+u(r)\exp (\mathrm{i}\lambda t+\mathrm{i}M\theta )+v^{\ast }(r)\exp (-\mathrm{i}{\lambda }^{\ast }t-\mathrm{i}M\theta )], where M is the azimuthal index carried by small perturbations u and v. The stationary solutions are linearly unstable if at least one pair of eigenvalue λ of the matrix operator \begin{array}{cc}{C}_{+}D& -D-C_{-}\end{array}, with C_{\pm }={\mathrm{d}}^2/\mathrm{d}{r}^2+r^{-1}\mathrm{d}/\mathrm{d}r-{(S\pm M)}^2/r^2-2g{\varphi }^2+\mu -V and D=-g{\varphi }^2, has an imaginary part. We have found that without the external potential, all exact VSs are subject to the azimuthal modulation instability. The exact VSs with n=2 may be stable in the presence of the harmonic potential, if the chemical potential is taken near the cutoff value {\mu }_{\mathrm{co}} (see Fig. 18). We stress that the mechanism supporting the stability in this case is different from the usual situation, when the soliton’s power is very small for the chemical potential approaching the cutoff value, which makes the NLSE an effectively linear equation [65].

4 Conclusions

In summary, we have examined the controlled dynamics of matter wave solitons in BECs with tunable interactions and external potential. We present the integrable conditions for the one and two-dimensional Gross-Pitaevskii equations, and obtain a family of exact analytical solutions. The dynamics of soliton, ring dark soliton, and vortex are investigated analytically and numerically. Our results show that the dynamics of the solitons can be effectively controlled by the Feshbach resonance which offers a good opportunity for the manipulation of atomic matter waves and nonlinear excitations in BECs. We hope that these results can stimulate studies of novel phenomena in nonlinear complex systems in general and further studies of BECs solitons in particular.

References

[1] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science, 1995, 269: 198

[2] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys., 1999, 71: 463

[3] P. G. Drazin and R. S. Johnson, Solitons: An Introduction, Cambridge: Cambridge University Press, 1989

[4] K. E. Strecker, G. B. Partridge, A. G. Truscoff, and R. G. Hulet, Nature (London), 2002, 417: 150

[5] L. Khaykovich, F. Schreck, G. Ferrari, T. Bourdel, J. Cubizolles, L. D. Carr, Y. Castin, and C. Salomon, Science, 2002, 296: 1290

[6] C. Becker, S. Stellmer, P. Soltan-Panahi, S. Dörscher, M. Baumert, Eva-Maria Richter, Jochen Kronjäger, K. Bongs, and K. Sengstock, Nature Phys., 2008, 4: 496

[7] B. Eiermann, Th. Anker, M. Albiez, M. Taglieber, P. Treutlein, K. P. Marzlin, and M. K. Oberthaler, Phys. Rev. Lett., 2004, 92: 230401

[8] D. L. Wang, X. H. Yan, and W. M. Liu, Phys. Rev. E, 2008, 78: 026606

[9] Z. W. Xie, W. P. Zhang, S. T. Chui, and W. M. Liu, Phys. Rev. A, 2004, 69: 053609

[10] Z. W. Xie, Z. X. Cao, E. I. Kats, and W. M. Liu, Phys. Rev. A, 2005, 71: 025601

[11] Z. D. Li and Q. Y. Li, Ann. Phys. (New York), 2007, 322: 1961

[12] Q. Y. Li and Z. D. Li, Opt. Commun., 2010, 283: 3361

[13] A. J. Moerdijk, B. J. Verhaar, and A. Axelsson, Phys. Rev. A, 1995, 51: 4852

[14] J. L. Roberts, N. R. Claussen, James P. Burke, Jr., C. H. Greene, E. A. Cornell, and C. E. Wieman, Phys. Rev. Lett., 1998, 81: 5109

[15] J. Stenger, S. Inouye, M. R. Andrews, H. -J. Miesner, D. M. Stamper-Kurn, and W. Ketterle, Phys. Rev. Lett., 1999, 82: 2422

[16] S. Inouye, M. R. Andrews, J. Stenger, H. -J. Miesner, D. M. Stamper-Kurn, and W. Ketterle, Nature, 1998, 392: 151

[17] S. L. Cornish, N. R. Claussen, J. L. Roberts, E. A. Cornell, and C. E. Wieman, Phys. Rev. Lett., 2000, 85: 1795

[18] E. A. Donley, N. R. Claussen, S. L. Cornish, J. L. Roberts, E. A. Cornell, and C. E. Wieman, Nature (London), 2001, 412: 295

[19] C. A. Regal and D. S. Jin, Phys. Rev. Lett., 2003, 90: 230404

[20] T. Volz, S. Dürr, S. Ernst, A. Marte, and G. Rempe, Phys. Rev. A, 2003, 68: 010702

[21] F. K. Abdullaev, A. M. Kamchatnov, V. V. Konotop, and V. A. Brazhnyi, Phys. Rev. Lett., 2003, 90: 230402

[22] V. M. Pérez-García, V. V. Konotop, and V. A. Brazhnyi, Phys. Rev. Lett., 2004, 92: 220403

[23] Z. X. Liang, Z. D. Zhang, and W. M. Liu, Phys. Rev. Lett., 2005, 94: 050402

[24] H. Saito and M. Ueda, Phys. Rev. Lett., 2003, 90: 040403

[25] D. E. Pelinovsky, P. G. Kevrekidis, and D. J. Frantzeskakis, Phys. Rev. Lett., 2003, 91: 240201

[26] G. D. Montesinos, H. Michinel, and V. M. Perez-Garcia, Phys. Rev. Lett., 2004, 92: 133901

[27] V. V. Konotop and P. Pacciani, Phys. Rev. Lett., 2005, 94: 240405

[28] G. Thalhammer, G. Barontini, L. De Sarlo, J. Catani, F. Minardi, and M. Inguscio, Phys. Rev. Lett., 2008, 100: 210402

[29] H. W. Xiong, S. J. Liu, W. P. Zhang, and M. S. Zhan, Phys. Rev. Lett., 2005, 95: 120401

[30] C. J. Myatt, E. A. Burt, R. W. Ghrist, E. A. Cornell, and C. E. Wieman, Phys. Rev. Lett., 1997, 78: 586

[31] D. S. Hall, M. R. Matthews, J. R. Ensher, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett., 1998, 81: 1539

[32] D. M. Stamper-Kurn, M. R. Andrews, A. P. Chikkatur, S. Inouye, H.-J. Miesner, J. Stenger, and W. Ketterle, Phys. Rev. Lett., 1998, 80: 2027

[33] T. S. Raju, P. K. Panigrahi, and K. Porsezian, Phys. Rev. A, 2005, 71: 035601

[34] G. Modugno, G. Ferrari, G. Roati, R. J. Brecha, A. Simoni, and M. Inguscio, Science, 2001, 294: 1320

[35] H. A. Cruz, V. A. Brazhnyi, V. V. Konotop, G. L. Alfimov, and M. Salerno, Phys. Rev. A, 2007, 76: 013603

[36] Q. Y. Li, Z. D. Li, S. F. Yao, L. Li, and G. S. Fu, Chin. Phys. B, 2010, 19: 080501

[37] V. M. Pérez-García and V. Vekslerchik, Phys. Rev. E, 2003, 67: 061804

[38] V. A. Brazhnyi and V. V. Konotop, Phys. Rev. E, 2005, 72: 026616

[39] J. K. Xue, G. Q. Li, A. X. Zhang, and P. Peng, Phys. Rev. E, 2008, 77: 016606

[40] L. Li, B. A. Malomed, D. Mihalache, and W. M. Liu, Phys. Rev. E, 2006, 73: 066610

[41] L. Salasnich and B. A. Malomed, Phys. Rev. A, 2006, 74: 053610

[42] P. Öhberg and L. Santos, J. Phys. B, 2001, 34: 4721

[43] S. K. Adhikari, Phys. Lett. A, 2005, 346: 179

[44] S. B. Papp, J. M. Pino, and C. E. Wieman, Phys. Rev. Lett., 2008, 101: 040402

[45] A. P. Sheppard and Y. S. Kivshar, Phys. Rev. E, 1997, 55: 4773

[46] D. J. Kaup and B. A. Malomed, Phys. Rev. A, 1993, 48: 599

[47] Th. Busch and J. R. Anglin, Phys. Rev. Lett., 2001, 87: 010401

[48] X. F. Zhang, Q. Yang, J. F. Zhang, X. Z. Chen, and W. M. Liu, Phys. Rev. A, 2008, 77: 023613

[49] C. Sulem and P. L. Sulem, The Nonlinear Schrödinger Equation, New York: Springer-Verlag, 1999

[50] B. Li, X. F. Zhang, Y. Q. Li, Y. Chen, and W. M. Liu, Phys. Rev. A, 2008, 78: 023608

[51] Y. Kagan, A. E. Muryshev, and G. V. Shlyapnikov, Phys. Rev. Lett., 1998, 81: 933

[52] S. V. Manakov, Sov. Phys. JETP, 1974, 38: 248

[53] R. Radhakrishnan and M. Lakshmanan, J. Phys. A, 1995, 28: 2683

[54] P. Öhberg and L. Santos, Phys. Rev. Lett., 2001, 86: 2918

[55] X. F. Zhang, X. H. Hu, X. X. Liu, and W. M. Liu, Phys. Rev. A, 2009, 79: 33630

[56] H. Pu and N. P. Bigelow, Phys. Rev. Lett., 1998, 80: 1130

[57] X. X. Liu, H. Pu, B. Xiong, W. M. Liu, and J. B. Gong, Phys. Rev. A, 2009, 79: 013423

[58] D. S. Wang, X. H. Hu, and W. M. Liu, Phys. Rev. A, 2010, 82: 023612

[59] K. Kasamatsu and M. Tsubota, Phys. Rev. A 74, 013617 (2006).

[60] X. H. Hu, X. F. Zhang, D. Zhao, H. G. Luo, and W. M. Liu, Phys. Rev. A, 2009, 79: 023619

[61] J. Belmonte-Beitia, V. M. Pérez-García, and V. Vekslerchik, Phys. Rev. Lett., 2007, 98: 064102

[62] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, New York: Dover, 1965

[63] D. S. Wang, X. H. Hu, J. P. Hu, and W. M. Liu, Phys. Rev. A, 2010, 81: 025604

[64] E. T. Whittaker and G. N. Watson, A Course in Modern Analysis, 4th Ed., Cambridge, UK: Cambridge University Press, 1990

[65] L. Wu, L. Li, J. F. Zhang, D. Mihalache, B. A. Malomed, and W. M. Liu, Phys. Rev. A, 2010, 81: 061805(R)

Rights & permissions

©Higher Education Press and Springer-Verlag Berlin Heidelberg 2010