Radially Symmetric Nonlinear States of Harmonically Trapped
Starting from the spectrum of the radially symmetric quantum harmonic oscillator in two dimensions, we create a large set of nonlinear solutions. The relevant three principal branches, with and radial nodes respectively, are systematically continued as a function of the chemical potential and their linear stability is analyzed in detail, in the absence as well as in the presence of topological charge , i.e., vorticity. It is found that for repulsive interatomic interactions only the ground state is linearly stable throughout the parameter range examined. Furthermore, this is true for topological charges or ; solutions with higher topological charge can be unstable even in that case. All higher excited states are found to be unstable in a wide parametric regime. However, for the focusing/attractive case the ground state with and can only be stable for a sufficiently low number of atoms. Once again, excited states are found to be generically unstable. For unstable profiles, the dynamical evolution of the corresponding branches is also followed to monitor the temporal development of the instability.
The study of trapped Bose-Einstein condensates (BECs) has had a high impact in recent years in a number of fields, including atomic, molecular, and optical physics, nuclear physics, condensed matter physics, chemical physics, applied mathematics, and nonlinear dynamics stringari (); pethick (); BECBOOK (). From the point of view of the latter, the topic of particular interest here is that at the mean field level, the inter-particle interaction is represented as a classical, but nonlinear self-action dalfovo (), leading to the now famous Gross-Pitaevskii (GP) equation as a celebrated model for BECs in appropriate settings. This has resulted in a large volume of studies of nonlinear excitations, including the prediction of excited states yukalov_add (), the experimental observation of dark dark (); dark1 (); dark2 (); dark3 (), bright expb1 (); expb2 () and gap gap () solitons in quasi-one-dimensional systems, as well as theoretical and experimental investigations of vortices, vortex lattices review1 (); review2 (), and ring solitons djf (); brand (); carr () in quasi-two-dimensional systems.
Apart from purely nonlinear dynamical techniques, such as the perturbation theory for solitons employed in Ref. djf (), other methods based on the corresponding linear Schrödinger problem may also be employed for the study of excited BECs. In particular, the underlying linear system for a harmonic external trapping potential is the quantum harmonic oscillator (QHO) stringari (); pethick (); dalfovo (), whose eigenstates are well-known since Schrödinger’s original treatment of the problem in 1926. Beginning with the linear equation, solutions can be numerically continued to encompass the presence of the nonlinear representation of the inter-particle interaction. Then, the QHO states can persist, bifurcate, or even disappear. This path does not seem to have been exploited in great detail in the literature. In one spatial dimension, it has been used in Ref. tristram (); jpb (), to illustrate the persistence tristram () and dynamical relevance jpb () of the nonlinear generalization of the QHO eigenstates. In higher dimensions, the work of Ref. tristram2 () illustrated the existence of solutions in a radially symmetric setting. Further progress has been hampered by the additional difficulties in (a) examining the linear stability and (b) converting the radial coordinate system to a Cartesian one to study evolution dynamics, including nonlinear stability. However, the mathematical tools for such an analysis exist, as we discuss in more detail below, and have to a considerable extent been used in Ref. carr (), especially in connection with ring-like structures with vorticity.
In this paper we present a systematic analysis of the existence, linear stability, and evolution dynamics of the states that exist in the spectrum of harmonically confined condensates from the linear limit, and persist in the nonlinear problem. Our analysis shows that in the case of repulsive interactions (defocusing nonlinearity), the only branch that is linearly stable consists of the ground state solution, i.e., a single hump with radial nodes. All higher excited states are linearly unstable and break up in ways that we elucidate below, if evolved dynamically in our system; typically this last stage involves also the loss of radial symmetry to unstable azimuthal perturbations. Here, we will treat explicitly the cases of and radial nodes; we have confirmed similar results for higher number of nodes.
On the other hand, in the case of attractive interactions (focusing nonlinearity), the system is subject to collapse in the absence of the potential. In fact, this constitutes the critical dimension for the underlying nonlinear problem sulem () beyond which collapse is possible. In this case, we observe that only the ground state branch with small values of the chemical potential, i.e., a small number of atoms, can be marginally stable. Once again, all excited states are unstable and the typical scenario here involves the manifestation of collapse-type catastrophic instabilities sulem (), as we elucidate through numerical simulations.
Our presentation is structured as follows: in section II, we provide the theoretical setup and numerical techniques. In section III, we present and discuss the relevant results. In section IV, we summarize our findings and present our conclusions. Finally, in the appendix we provide relevant details regarding our numerical methods.
Ii Setup and Numerical Methods
We will use as our theoretical model the well-known GP equation in a two-dimensional (2D) setting. As is well-known, the “effective” 2D GP equation applies to situations where the condensate has a nearly planar, i.e., “pancake” shape (see, e.g., Ref. Beersheva () and references therein.) We express the equation in harmonic oscillator units rupr () in the form:
where is the 2D wave function (the subscript denotes partial derivative). The wave function is the condensate order parameter, and has a straightforward physical interpretation: is the local condensate density. The external potential assumes the typical harmonic form, with being the effective frequency of the parabolic trap; the latter can be expressed as , where are the confinement frequencies in the radial and axial directions, respectively. It is assumed that for the pancake geometry considered herein; in particular, we choose for our computations. Finally, is the normalized coefficient of the nonlinear term, which is fixed to for attractive interactions and to for repulsive interactions. Accordingly, the squared norm is
where is the number of atoms and is the usual nonlinear coefficient in the GP equation in harmonic oscillator units, renormalized to two dimensions appropriately carr2005c ().
In order to numerically identify stationary nonlinear solutions of Eq. (1), we use the standing wave ansatz , where we assume the density of the solution is radially symmetric with chemical potential and topological charge (vorticity) . This results in the steady-state equation:
Equation (3) exhibits infinite branches of nonlinear bound states, each branch stemming from the corresponding mode of the underlying linear problem. These branches are constructed and followed using the method of pseudo-arclength continuation arclength1 (); arclength2 (), in order to obtain subsequent points along a branch, once a point on it has been identified. The first point on each branch is found using a bound state of the underlying linear equation
The linear state corresponds to a solution of Eq. (3) with parameters , frequency , vorticity , number of radial nodes and normalization tending to zero; in that limit, the nonlinearity becomes negligible and the linear solutions are the well-known Gauss-Laguerre modes with . Using such an initial guess, a non-trivial solution is found through a fixed-point iterative scheme for a slightly perturbed value of the chemical potential. This is done on a grid of Chebyshev points suited to the radial problem, following the approach of Ref. trefethen () for the Laplacian part of the equation, as is explained in detail in the appendix. The pseudo-arclength method, used to trace subsequent solutions, works via the introduction of a pseudo-arclength parameter and an additional equation such that where is a solution of Eq. (3). We used for :
Lastly, the new expanded system of equations is solved using a predictor-corrector method.
The linear stability of the solutions is analyzed by using the following ansatz for the perturbation:
where the asterisk stands for the complex conjugation and is the polar angle. One solves the resulting linearized equations for the perturbation eigenmodes and eigenvalues associated with them. The key observation that allows one to carry this task through is that the subspace of the different azimuthal perturbation eigenmodes, i.e., subspaces of different , decouple pego (); carr (), leading each eigenmode to be coupled only with its complex conjugate. This in turn allows the examination of the stability of the 2D problem in the form of a denumerably infinite set of one-dimensional radial eigenvalue problems that are solved on the same grid of Chebyshev points as the original existence problem of Eq. (3).
In order to examine the dynamics of the cases found to be unstable, to confirm the validity of the cases identified as stable, and to test for nonlinear instabilities, the 1D radial solutions of Eq. (3) were taken as an initial condition for a code which solved Eq. (1) with a Chebyshev spectral radial-polar method in space and a fourth-order integrator in time. The Chebyshev spectral radial-polar method is the most natural choice for our setting, as it a) avoids the conversion of the radial solution into an interpolated Cartesian grid, and, more importantly, b) avoids spurious effects associated with a mismatch between the symmetry of the solution and that of the grid. For example, we have observed that the use of a Cartesian grid artificially enhances the excitation of modes that have a similar symmetry as the grid, and in particular the mode. The results of all of the above numerical techniques, i.e., existence, linear stability, and dynamical evolution, are reported in what follows.
Our steady state and stability results are summarized in Fig. 1. The panel shows the continuation from the linear limit of the relevant states. The nonlinearity leads to a decrease of the chemical potential with increasing power for the focusing case and to a corresponding increase for the defocusing case. Notice that some of the branches presented here have also appeared in the earlier work of Ref. tristram2 (). However, the important new ingredient of the present work is the detailed examination, both through linear stability analysis and through dynamical evolution, of the stability of these solutions. The latter is encapsulated straightforwardly in the symbols associated with each branch. The plus symbols denote branches which are stable, while the diamond symbols correspond to unstable cases. We observe that the only truly stable solution corresponds to the ground state with of the defocusing problem, which for large number of atoms can be well approximated by the Thomas-Fermi state stringari (); pethick (). Furthermore, this is true only for topological charges and ; for higher topological charges , there are regions of stability for , as was originally shown in Refs. pu () and vortex_stability_add (). All higher excited states with one, two or more nodes are directly found to be unstable, even close to the corresponding linear limit of these states. In the focusing case, the situation is even more unstable due to the catastrophic effects of self-focusing and wave collapse sulem (); doddRJ1996 (). In particular, even the ground state with may be unstable due to collapse (represented by the mode with in this case), although the instability growth rate may be very small, as is the case for , see e.g., Fig. 2. Excited states are always unstable for the focusing nonlinearity also; in fact, they are more strongly so than in the defocusing case, again due to the presence of the mode.
Having offered an outline of the stability properties of our solutions in Fig. 1, we now examine the subject in detail in the following results. In Figs. 2–4 we depict the real part of the primary () eigenvalues for , , and , respectively and for the states with and in each case, i.e., the ground state and first two excited states in the top, middle and bottom panel of each figure. The stability results showcase the presence of at least one unstable eigenvalue. In fact, there is only one instability eigenvalue for the ground state of the attractive case, and a larger number of such eigenvalues for excited states. The presence of such eigenvalues, whose larger magnitude is tantamount to more rapid dynamical development of the instability, indicates that the dynamical evolution of such stationary states, when small perturbations are added to them, will manifest the presence of such unstable eigenmodes through the deformation and likely destruction of the initial structure.
This is shown in detail for in Figs. 5–7. In the top row of panels in Fig. 5 we depict the radial profiles of the solution for and (left panel) and and (right panel) — for all other figures we use when . The middle row of panels in the figure depicts their corresponding linear stability spectra in the complex plane . The unstable eigenvalues for different values of are given different symbols as described in the table in Table 1. Eigenvalues with or for are plotted with a small dark dot.
Finally, the bottom row of panels in Fig. 5 depicts the corresponding time evolution of the chosen profile after a random perturbation of amplitude was added to the steady state profile at time . For the profile considered in Fig. 5, namely the ground state with , it is clear that the solution for the repulsive case () is stable since it corresponds to the Thomas-Fermi ground state. On the other hand, the solution for the attractive case () is weakly unstable, due to the mode, as discussed above. The eigenvalues for this mode are depicted by the circles in the middle-left panel. This instability is due to the well-known collapse of the solution for the attractive case, where the solution is seen to tend to a thin spike carrying all the mass (see the time evolution in the bottom-left panels).
In Figs. 6 and 7 we present the equivalent results for the first and second excited states ( and respectively) and for the same case. In these cases the instability manifests itself with the presence of a richer scenario of unstable eigenvalues. Let us explain in detail the first excited state in Fig. 6. The first excited state in the attractive case (left panels) is still prone to collapse as evidenced by the quartet of eigenvalues depicted by the circles in the middle panel. This quartet is responsible for the collapse of the central spike of the solution as seen in the time series evolution (bottom panels). Furthermore, the and modes ( and symbols in the figure) are the most unstable ones and are responsible for the azimuthal modulational instability of the first ring of the solution as evident in the time evolution (bottom panels).
On the other hand, for the repulsive case (right panels) of the first excited state with it is clear that the eigenvalue is stable as there is no collapse in the repulsive case. The dynamic evolution of the instability in this case leads to the competition between the unstable eigenvalues , and (in order of strength) that is seen to be dominated by the most unstable mode , as can be evidenced by the three humps (surrounding the central peak) displayed at .
Finally, in Fig. 7 we depict similar results for the second excited state () for . As can be seen from the figures, the more excited the state, the richer the (in)stability spectra. It is worth noticing again that the attractive case is, as before, prone to collapse due to a strong unstable mode while, naturally, the repulsive case lacks this collapsing mode. The richer set of eigenvalues for higher excited states is easy to interpret since higher excited states possess more radial nodal rings.
In general, for the repulsive cases where collapse is absent, we observe that in all cases the solution is subject to azimuthal modulations that produce coherent structures reminiscent of the “azimuthons” of Ref. azimuthons (). On the other hand, in the attractive cases collapse is ubiquitous, due to the mode with ; in addition, azimuthal modulations emerge.
The above discussion describes in detail the stability and dynamics for the non-topologically charged () solutions. Let us now describe the case when the solutions have an intrinsic topological charge, namely . In Figs. 8–10 we display the results for the ground state, first and second excited states ( and radial nodes) with unit topological charge (or vorticity), namely . In fact, Fig. 8 corresponds to the singly charged vortex solution at the center of the harmonic trap. As it is well known, this solution is stable in the repulsive case (see right panels) and unstable in the attractive case (see left panels). It is interesting to note that the instability of the vortex solution in the attractive case is not driven by collapse ( eigenvalue) since the vorticity tends to push the mass away from the center () of the trap. This is a general feature of the cases with , where collapse appears to arise in the rings of the cloud, rather than at its center. Finally, in Figs. 11-13 we display similar results for the doubly charged case .
It is worth mentioning that the unstable modes for vorticity-carrying solutions () tend to rotate as they are growing. This effect can be clearly seen in the three-dimensional iso-density contour plot depicted in the top panel of Fig. 14 corresponding to the first excited state with and . All modes tend to stop rotating after the spikes created reach a certain height, as can be seen clearly in the bottom panel in Fig. 14. This figure shows the ground state with and . Interestingly, a single solution with multiple radial nodes (excited states) is able to pick out more than one growing mode since each ring can be affected by a different -mode. This effect is seen in the top-left panels in Fig. 12 where the first excited state for and is seen to develop, at earlier times, the mode in the inner ring while, at later times, the outer ring develops the mode.
In this paper, we have revisited the topic of nonlinear continuation of linear, radially-dependent Laguerre-Gauss states of the two-dimensional quantum harmonic oscillator model in the presence of inter-particle interaction-induced nonlinearity. We have systematically constructed such solutions starting from the linear limit and, more importantly, we have detailed their linear and nonlinear stability properties. This was accomplished by careful examination of the corresponding eigenvalue problem or, more appropriately, the (infinite) one-parameter family of eigenvalue problems, in radial coordinates. We have also provided a full numerical time evolution of the model on a radial-polar grid. We have principally observed that the ground state of the attractive case is unstable due to collapse, although the growth rate of the instability may be weak. For the repulsive case the relevant state is stable for topological charges and ; for higher charges the stability depends on the atom number. Excited states have been found to be generically unstable in both attractive and repulsive cases; the development of the instabilities produces collapse in the former, while it results in the formation of azimuthally modulated states in the latter.
It would certainly be of interest to extend the present techniques to the full 3D problem, again considering radial states and their continuation from the linear limit, as well as their linear stability. However, in the latter setting direct numerical computations are significantly more intensive. Such studies are outside of the scope of the present work and will be considered in a future publication.
P.G.K. gratefully acknowledges the support of NSF-DMS-0204585 and NSF-CAREER. P.G.K. and R.C.G. acknowledge the support of NSF-DMS-0505663, and L.D.C. acknowledges support under NSF grant PHY-0547845 as part of the NSF CAREER program. L.D.C. acknowledges useful discussions with Charles Clark and Mark Edwards.
Appendix: Spectral Methods
When dealing with the Laplacian in polar coordinates, the origin presents a significant challenge due to division by zero:
This problem has been faced with in the past tristram (); tristram2 (); jpb (); carr (), with respect to the GP equation, but the methods used have so far limited the scope of such investigations. However, spectral methods can be used to avoid the singularity to handle the Laplacian in studies of the GP equation.
The particular application of spectral methods we used in our calculations is described by Trefethen in Ref. trefethen (), but is restated here for completeness. This method avoids the problem with the origin by avoiding the origin altogether through the use of an even number of Chebyshev nodes. In order to check for the accuracy of utilizing spectral methods for the purposes of modeling the GP equation in polar coordinates, we recreated the results of Fig. 1b in Ref. pu (). As can be seen from Fig. 15, our spectral method code has generated a plot in close agreement with the original plot.
Spectrally discretizing the polar domain is achieved by discretizing the radial direction using polynomial interpolation, in particular the Chebyshev nodes:
meaning the domain must be normalized in order for the solution to be properly contained within , while the angular direction is discretized using Fourier (or trigonometric) interpolation:
These choices for discretization are based upon the periodicity of the angular direction and the lack thereof in the radial direction (which causes the Gibbs phenomenon to occur if Fourier interpolation is used). Additionally, the Chebyshev polynomials
Implementation is simply done by using the nodes described above along with their corresponding derivative matrices. For the Chebyshev nodes, the Chebyshev differentiation matrix is given by:
while the Fourier differentiation matrix (for an even number of Fourier nodes) is given by:
It should be noted that the differentiation matrix for an odd number of Fourier nodes differs from the one above due to a correction which must be made in the derivation of the differentiation matrix for an even number of Fourier nodes trefethen (), but we have chosen to restrict ourselves to using even numbers of Fourier nodes.
Typically when using polar coordinates, the domain of the radial direction is restricted to since any function on the unit disc must inherently satisfy the symmetry condition:
and since we have chosen to use an even number of Fourier nodes, our grid must also satisfy this condition. As a result, solution values would be recorded and evaluated twice using this grid. Initially, this may appear to be a significant disadvantage for this discretization, but an elegant simplification can be arrived at using a symmetry property of the Chebyshev spectral differentiation matrix: .
Using both symmetries, the derivative of in the radial direction can be shown to be exactly the same for both occurrences of the node:
In addition, the derivative calculation can be restricted to using only the positive -axis by using (assuming ):
The same simplification can be inductively shown to work for higher derivatives since they are calculated by multiple applications of the derivative matrix
and Eq. (Appendix: Spectral Methods) has already established the symmetry of the first derivative, thus by induction, all higher derivatives can be calculated using only the positive -axis. Consequently, any simulations based upon a Chebyshev-Fourier polar grid can be written to strictly use the positive -axis, resolving the computational time and storage problems, but still benefiting from the advantages of spectral methods. (See Chapter 11 of Ref. trefethen () for an implementation of this method in Matlab.)
Spectral methods can also be used to evaluate integrals with greater accuracy. In particular, this was used for calculating the power of the numerically determined steady state solutions. To further explore how spectral methods can be used for integration, take a generic integral
for some sufficiently smooth function . This integral can then be restated in the form of an initial value ODE
and discretized using Chebyshev nodes, thus replacing the derivative by the Chebyshev spectral differentiation matrix
where here is the discrete version of the integrand above, , and the boundary condition, , is imposed by stripping off the last row and column of to obtain the matrix . The value of the integral at every Chebyshev node is then easily found by inverting the derivative matrix to get . However, the only value we are interested in is (i.e. the integral over the entire domain) and this can easily be found by using only the first row of , call it : .
- (1) L.P. Pitaevskii and S. Stringari, Bose-Einstein Condensation, Oxford University Press (Oxford, 2003).
- (2) C.J. Pethick and H. Smith, Bose-Einstein condensation in dilute gases, Cambridge University Press (Cambridge, 2002).
- (3) P.G. Kevrekidis, D.J. Frantzeskakis, and R. Carretero-González (eds). Emergent Nonlinear Phenomena in Bose-Einstein Condensates: Theory and Experiment. Springer Series on Atomic, Optical, and Plasma Physics, Vol. 45, 2008.
- (4) F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999).
- (5) V.I. Yukalov, E.P. Yukalova, and V.S. Bagnato, Phys. Rev. A 56, 4845 (1997); P.W. Courteille, V.S. Bagnato, and V.I. Yukalov, Laser Phys. 11, 659 (2001).
- (6) S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Sengstock, A. Sanpera, G.V. Shlyapnikov, and M. Lewenstein, Phys. Rev. Lett. 83, 5198 (1999);
- (7) J. Denschlag, 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 287, 97 (2000);
- (8) B.P. Anderson, P.C. Haljan, C.A. Regal, D.L. Feder, L.A. Collins, C.W. Clark, and E.A. Cornell, Phys. Rev. Lett. 86, 2926 (2001);
- (9) Z. Dutton, M. Budde, Ch. Slowe, and L.V. Hau, Science 293, 663 (2001).
- (10) K. E. Strecker, G. B. Partridge, A. G. Truscott, and R. G. Hulet, Nature 417, 150 (2002).
- (11) L. Khaykovich, F. Schreck, G. Ferrari, T. Bourdel, J. Cubizolles, L. D. Carr, Y. Castin, and C. Salomon, Science 296, 1290 (2002).
- (12) B. Eiermann, Th. Anker, M. Albiez, M. Taglieber, P. Treutlein, K.-P. Marzlin, and M. K. Oberthaler, Phys. Rev. Lett. 92, 230401 (2004).
- (13) A.L. Fetter and A.A. Svidzinsky, J. Phys.: Cond. Matt. 13, R135 (2001).
- (14) P.G. Kevrekidis, R. Carretero-González, D.J. Frantzeskakis and I.G. Kevrekidis, Mod. Phys. Lett. B 18, 1481 (2004).
- (15) G. Theocharis, D.J. Frantzeskakis, P.G. Kevrekidis, B.A. Malomed and Yu.S. Kivshar, Phys. Rev. Lett. 90, 120403 (2003).
- (16) N.S. Ginsberg, J. Brand and L.V. Hau, Phys. Rev. Lett. 94, 040403 (2005).
- (17) L.D. Carr and C.W. Clark, Phys. Rev. Lett. 97, 010403 (2006); Phys. Rev. A 74, 043613 (2006).
- (18) Yu.S. Kivshar, T.J. Alexander and S.K. Turitsyn, Phys. Lett. A 278, 225 (2001).
- (19) P.G. Kevrekidis, V.V. Konotop, A. Rodrigues and D.J. Frantzeskakis, J. Phys. B: At. Mol. Opt. Phys. 38, 1173 (2005).
- (20) Yu.S. Kivshar and T.J. Alexander, cond-mat/9905048;
- (21) C. Sulem and P.L. Sulem, The Nonlinear Schrödinger Equation, Springer-Verlag (New York, 1999).
- (22) Y.B. Band, I. Towers, and B.A. Malomed, Phys. Rev. A 67, 023602 (2003).
- (23) P.A. Ruprecht, M.J. Holland, K. Burnett and M. Edwards, Phys. Rev. A 51, 4704 (1995).
- (24) L. D. Carr, M. J. Holland, and B. A. Malomed, J. Phys. B: At. Mol. Opt. 38, 3217 (2005).
- (25) E. L. Allgower and K. Georg, Numerical Path Following. In P. G. Ciarlet and J. L. Lions, Handbook of Numerical Analysis, Volume 5 3-207 (North-Holland (1997)).
- (26) H.B. Keller, Numerical Solution of Bifurcation and Nonlinear Eigenvalue Problems. In P.H. Rabinowitz, Applications of Bifurcation Theory 359-384 (Academic Press, New York, 1977).
- (27) L.N. Trefethen, Spectral Methods in Matlab, SIAM (Philadelphia, 2000).
- (28) R.L. Pego and H. Warchall, J. Nonl. Sci. 12 (2002) 347–394.
- (29) H. Pu, C.K. Law, J.H. Eberly, and N.P. Bigelow, Phys. Rev. A 59, 1533 (1999).
- (30) D. Mihalache, D. Mazilu, B.A. Malomed, and F. Lederer, Phys. Rev. A 73, 043615 (2006); B.A. Malomed, F. Lederer, D. Mazilu, and D. Mihalache, Phys. Lett. A 361, 336 (2007); J.A.M. Huhtamäki, M. Möttönen, and S.M.M. Virtanen, Phys. Rev. A 74, 063619 (2007); E. Lundh and H.M. Nilsen, Phys. Rev. A 74, 063620 (2007).
- (31) R. J. Dodd, M. Edwards, C. J. Williams, C. W. Clark, M. J. Holland, P. A. Ruprecht, and K. Burnett, Phys. Rev. A 54, 661 (1996).
- (32) A.S. Desyatnikov, A.A. Sukhorukov and Yu.S. Kivshar, Phys. Rev. Lett. 95, 203904 (2005).
- (33) R. Peyret, Spectral Methods for Incompressible Viscous Flow, Springer-Verlag (New York, 2002).