Emergence and Stability of Vortex Clusters in BoseEinstein Condensates:
a Bifurcation Approach near the Linear Limit
Abstract
We study the existence and stability properties of clusters of alternating charge vortices in BoseEinstein condensates. It is illustrated that such states emerge from cascades of symmetrybreaking bifurcations that can be analytically tracked near the linear limit of the system via weakly nonlinear fewmode expansions. We present the resulting states that emerge near the first few eigenvalues of the linear limit, and illustrate how the nature of the bifurcations can be used to understand their stability. Rectilinear, polygonal and diagonal vortex clusters are only some of the obtained states while mixed states, consisting of dark solitons and vortex clusters, are identified as well.
keywords:
BoseEinstein condensates, Vortices, Dark solitons, BifurcationsPacs:
03.75.Lm, 67.90.+z, 34.50.Cx, , , , and
1 Introduction
Vortices are among the most striking characteristics of nonlinear field theories in higherdimensional settings [1]. They constitute one of the remarkable features of superfluids, while playing also a key role in critical current densities and resistances of typeII superconductors through their transport properties, and are associated with quantum turbulence in superfluid helium [2]. Vortices appear also in a wide variety of fields, ranging from fluid dynamics [3] to atomic physics [4] and optical physics [5].
A pristine setting for the study of vortices at the mesoscale has emerged after the realization of atomic BoseEinstein condensates (BECs). In this context, socalled matterwave vortices were experimentally observed therein [6], by using a phaseimprinting method between two hyperfine spin states of a Rb BEC [7]. This achievement subsequently triggered extensive studies concerning vortex formation, dynamics and interactions. For instance, stirring the BECs [8] above a certain critical angular speed [9, 10, 11, 12] led to the production of few vortices [12], and even of robust vortex lattices [13]. Vortices can also be formed in experiments by means of other techniques, such as by dragging obstacles through the BEC [14] or by the nonlinear interference of different condensate fragments [15]. Not only unitcharged, but also highercharged vortices were observed [16] and their dynamical instabilities have been analyzed.
While a considerable volume of work has been dedicated to individual vortices and to vortex lattices, arguably, vortex clusters consisting of only a few vortices have attracted less interest. The latter theme has become a focal point recently, through the experiments involving twovortex states (alias vortex dipoles) [17, 18], as well as threevortex states [19]. In Ref. [17], vortex dipoles were produced by dragging a localized light beam with appropriate speed through the BEC, while in Ref. [18] they were distilled through the KibbleZurek mechanism [20], previously proposed and realized for vortices in Ref. [21]. For the nonlinear dynamics of the vortices in the dipoles of Ref. [18], see also the very recent analysis of Ref. [22]. In Ref. [19], different types of threevortex configurations were produced by applying an external quadrupolar magnetic field on the BEC. The principal ones among them were an aligned vortex “tripole” with a vortex of one topological charge straddling two other oppositely charged vortices, and an equilateral triangle of three same charge vortices. On the theoretical side, few vortex states have been considered also in a number of works. It was shown, in particular, that vortex dipoles (consisting of a pair of vortices with opposite circulation) are fairly robust in BECs [23]. More elaborate states, such as dipoles, tripoles and quadrupoles, were considered in Refs. [24, 25]. Dynamics of such few vortex states in the weaklyinteracting limit were performed in Ref. [26], while the recent work of Ref. [27] connected the vortex dipoles to the instability of dark soliton stripes; see also the important earlier work of Ref. [28].
In this work, we present a unifying analysis of the existence and stability of vortex clusters (consisting of alternating charge vortices), by corroborating theoretical investigations and numerical computations. Our study in Section 2 will be based on the lowdensity limit of nearlinear excitations, where we will illustrate how they emerge (bifurcate through symmetrybreaking bifurcations) from states of the twodimensional (2D) quantum harmonic oscillator. Then, such theoretically identified states will be continued via numerical computations in Section 3 to the strongly nonlinear regime. From these continuations, we will be able to infer numerous previously undiscovered vortex cluster states, and to elucidate their stability properties (as well as compare to the theoretical predictions). Lastly, in Section 4, we will summarize our findings and present some directions for future work.
2 Model And Theoretical Analysis
We consider a quasi2D (alias “diskshaped”) condensate confined in a highly anisotropic trap with frequencies and along the transverse and inplane directions, respectively. In the case and (where is the chemical potential), and for sufficiently low temperatures, the inplane part of the macroscopic BEC wave function obeys the following dimensional GrossPitaevskii equation (GPE) (see, e.g., Ref. [4]):
(1) 
where is the inplane Laplacian, while the potential is given by (where is the atomic mass). The effective 2D nonlinearity strength is given by , with , and denoting, respectively, the threedimensional (3D) interaction strength, the wave scattering length, and the transverse harmonic oscillator length. Equation (1) can be expressed in the following dimensionless form,
(2) 
where the density , length, time and energy are respectively measured in units of , , and . Finally, the harmonic potential is now given by , with . From here on, all equations will be presented in dimensionless units for simplicity.
Below, we will analyze the existence and linear stability of the nonlinear modes of Eq. (2). Notice that numerically the relevant nonlinear states will be identified as a function of the chemical potential by means of a fixed point (Newton iteration) scheme over a rectangular twodimensional grid with suitably small spacing. We will also explore the linear (spectral) stability of the obtained states by means of the Bogoliubovde Gennes (BdG) analysis. The latter involves the derivation of the BdG equations, which stem from a linearization of the GPE (2) around the stationary solution via the ansatz
(3) 
where denotes complex conjugate. The solution of the ensuing BdG eigenvalue problem yields the eigenfunctions and eigenfrequencies . Due to the Hamiltonian nature of the system, if is an eigenfrequency of the Bogoliubov spectrum, so are , and . Notice that a linearly stable configuration is tantamount to , i.e., all eigenfrequencies being real. It is important to mention that in what follows we only resolve the relevant eigenvalues up to due to computational domain constraints and therefore all eigenvalues smaller than will be omitted in the figures.
Within the BdG analysis, a relevant quantity to consider is the norm energy product of a normal mode with eigenfrequency , namely,
(4) 
The sign of this quantity, known as Krein sign [29], is a topological property of each eigenmode. In particular, if this sign is negative and such a mode becomes resonant with a mode with positive Krein signature then, typically, complex frequencies appear in the excitation spectrum, i.e., a dynamical instability arises [29]. We refer to this as an oscillatory instability. Furthermore, dynamical instabilities may arise due to a real mode eigenfrequency becoming imaginary. This typically coincides with a bifurcation of a new state.
In the context of Eq. (2), it is useful to consider the lowdensity (linear) limit. There, eigenstates of the 2D quantum harmonic oscillator arise in the form (as well as linear combinations thereof), where and , quantify the order (and number of nodal lines) in each direction. This produces a linear limit whose first excited state has and linear eigenstates and . Notice that one of their interesting linear combinations is , which creates the singlecharged vortex even at this linear limit; this state exists for all higher values of and is shown in Fig. 1, but we will not be concerned with it further herein, as our focus will be on clusters of vortices.


Each of the above mentioned linear eigenstates, and , represents a “stripe” i.e., a state with a nodal line. These can be continued for higher chemical potentials . As increases, this state develops into a onedimensional (1D) dark soliton stripe, which is an exact analytical solution of Eq. (2) in the absence of the trap [30]. However, it is wellknown that such a state is dynamically unstable towards decay into vortex structures [31, 32, 33]. This decay can be understood from a symmetrybreaking bifurcation point of view [27, 34]. In particular, it is possible to consider a twomode (Galerkintype) expansion, similar to the one used in the literature of doublewell potentials (see e.g. Ref. [35]) in the form:
(5) 
where , are complex timedependent prefactors, while , and . The resulting equations and analysis are formally equivalent to the ones derived in Ref. [35] (see Eqs. (4)(5) therein), with appropriate modifications of the inner products, but also with a fundamental difference. In the 1D doublewell setting, only symmetrybreaking bifurcations of asymmetric real solutions are predicted (the socalled states that have recently been experimentally observed in Ref. [36]). The richer 2D case enables bifurcations even when the relative phase between the complex order parameters and is . In particular, such bifurcations are generically predicted at an atom number:
(6) 
where , are the linear state eigenvalues corresponding to and , while and ; the critical chemical potential is given by .
3 Numerical Results and Comparison with Theory
The above mentioned twomode theory provides explicit predictions for the bifurcation not only of the vortex dipole (vd) state when , but also for an aligned vortex tripole (3v) state (cf. the experimental observations of Ref. [19]) for , for an aligned vortex quadrupole (4v) state with , etc. In fact, there is an entire cascade of such bifurcations, as increases, which occur progressively at , , , and for etc., respectively. Notice that the number of vortices of the resulting cluster is evident by the number of intersections of the single nodal line of with the perpendicular nodal lines of , as well as the relative phase of their complex prefactors at these intersections. Also, it is evident that that the sign changing of at these intersections leads to an alternation of the ensuing vortex charges. Importantly, general bifurcation theory can be used to identify the stability characteristics of the resulting states. In particular, since the stripe is dynamically stable as it emerges from the linear limit, the vortex dipole state that arises from it upon the first symmetrybreaking “event” () should inherit this stability. However, now, once the stripe has become unstable, all higher bifurcations with will necessarily result into dynamically unstable states.


Numerical results on the symmetrybreaking bifurcations resulting in the emergence of vortex cluster states from the first excited state (single dark soliton stripe) are summarized in Fig. 1 for . The emergence of () states can be observed to occur respectively at while the corresponding theoretical predictions are . Clearly, the twomode approach captures the fundamental phenomenology, although a slight progressive degradation of the agreement on the critical point arises due to the increasing departure from the linear limit. Additionally, the linear stability properties of the resulting states directly reflect the theoretical stability expectations discussed above; see also Ref. [34].
Importantly, this approach is not restricted to the first excited state. The advantage of the bifurcation method and of the wealth of states that can be derived from it is unveiled, e.g., when considering the next set of excited states, namely the combinations with , i.e., the degenerate states , and . In this setting, already a vortex quadrupole can be formed at the linear limit as [37]. Interestingly, so can a doublycharged vortex through . However, these states do not present symmetry breaking bifurcations and, therefore, are not considered further in what follows. Focusing on the states that do, some prototypical examples are (i) the solitonic state consisting of two stripes (i.e., a twodarksoliton state), i.e., the state, (ii) an Xshaped dark soliton cross emerging from , as well as (iii) a ring dark soliton state [38] (see also Refs. [39, 40]) arising from . The state is one for which the generalization of the phenomenology of Fig. 1 is most straightforward as with and states with (2 lines of vortices each) are formed with . For example, for and , such bifurcations are predicted at , and (see thick black vertical dashed lines in Fig. 2), respectively, leading to and twoline vortex clusters. On the other hand, the Xshaped dark soliton pair is subject to similar symmetrybreaking bifurcations/destabilizations due to , , etc. These bifurcations, occurring in turn at , , and (see thin orange [gray in printed version] vertical dashed lines in Fig. 2), lead to a diagonal state with 6 vortices, a doublycharged vortex in the center together with a fourvortex quadrupole around it, and an 8vortex cluster of neardiagonal vortices. Finally, the ring dark soliton gets mixed with states of the form , where , again for . Remarkably, these symmetrybreaking events, which can be predicted to occur, e.g., at and (see, respectively, orange [gray in printed version] and blue [dark in printed version] vertical dashed lines in Fig. 3), for and , give rise to polygonal vortex configurations with the vortices now placed on the periphery of the circle. This way, vortex hexagons, octagons, decagons, etc. can be systematically constructed at will.




In Figs. 2 and 3, we explore numerically the above symmetrybreaking events for the second excited states and show besides the bifurcation diagram also typical profiles of the waveforms, as they result from varying the chemical potential. Specifically, Fig. 2 depicts the bifurcations and profiles for the states arising from the two dark soliton stripes (middle column of panels) and the dark soliton cross (right column of panels). From the twostripe soliton bifurcates a sixvortex (6v2) () state and an eightvortex (8v2) state () in the considered regime. From the Xshaped dark soliton cross bifurcates a diagonal six (6x) and eightvortex (8x) configurations, as well as a state with a vortex quadrupole surrounding a central doublycharged vortex (5x). We also mention in passing that interesting additional bifurcation events (collisions and disappearances into “blue sky” bifurcations) arise, e.g., between the sixvortex cluster deriving from the Xshaped dark soliton cross and that deriving from the twostripe soliton (see blue circle denoted by A in Fig. 2). On the other hand, in Fig. 3 we depict the states that bifurcate from the vortex ring: vortex hexagons (6r) and octagons (8r). In this figure we also depict the vortex quadrupole state (bottom right panel) which does not bifurcate from the soliton ring. In all the cases, good quantitative agreement is found on the critical points predicted by the theory and those observed numerically.
Figures 4 and 5 complement the above existence picture with a systematic linear stability analysis of each of the corresponding states for Figs. 2 and 3 respectively. The twostripe soliton is dynamically stable near the linear limit, a stability inherited by its first symmetrybreaking offspring, the sixvortex () state; while the aligned eightvortex case () is unstable from its existence onset (see left column of panels in Fig. 4). Similarly, the Xshaped dark soliton cross state also bears imaginary eigenfrequencies for all values of and hence its derivative states inherit its dynamical instability (see right column of panels in Fig. 4). On the other hand, the ring dark soliton is, as was also found earlier [37, 38, 39, 40], always unstable, hence the polygonal vortices that derive from it inherit this instability (see Fig. 5). However, it should be noted that the instability of such states weakens as the chemical potential increases. Note that, the vortex quadrupole is stable in the linear limit. Since no state bifurcates from the vortex quadrupole this stability is generally preserved (apart from a small window of oscillatory instabilities).
A general comment on the observed BdG spectra is due here. We can see that the spectra bear a large number of positive Krein sign modes (see blue (black) points is all BdG spectra) which appear to be asymptoting to appropriate corresponding values in the large chemical potential limit. In addition, there is a number of negative Krein sign eigenfrequencies (see orange (dark gray) points is all BdG spectra) which may lead to oscillatory (upon collision with positive sign ones) or exponential instabilities. Generally, we can comment that the positive Krein sign eigenmodes correspond to the ground state “background” on which a particular solitonic or vortex (or mixed) state exists. Furthermore, these eigenfrequencies have a welldefined limit when is large as discussed e.g. in Ref. [41] (see also references therein). On the other hand, the negative Krein sign modes reflect the excited state nature of the considered soliton or vortex (or mixed) states and are the ones which bear the potential for dynamical instabilities. While for solitonic states, we do not have a precise count on the number of the latter eigenmodes, in the case of multivortex clusters consisting of individual vortices an upper bound on the maximal order of such (negative Krein) eigenstates can be given by the number of vortices in the configuration. This, in turn, also gives an upper bound on the number of potentially observable unstable eigenmodes.
Importantly, the present approach can be generalized to higher excited states. As merely a small sample of further exotic configurations that can emerge (some of which can even be structurally robust), we mention a 9vortex cluster (a square grid of vortex “particles”), which emerges from the linear limit as and can, therefore, be regarded as a higher excited analog to the single vortex and the vortex quadrupole states and is expected to be stable in the vicinity of that limit (at least with respect to purely imaginary [i.e., nonoscillatory] instabilities, see below). On the other hand, there are bifurcations from that state including the bifurcation of a ring with a vortex in its center (which was again discussed e.g., in Ref. [40]). Another state that exists and should be robust near the linear limit is a threesolitonstripe. In fact, from such a state bifurcations again emerge due to the mixture with states such as with , but also with with . We focus briefly on the latter, which is theoretically predicted to occur at (see vertical dashed line in Fig. 6), because it gives rise to yet another novel type of state, namely a mixed state between vortices and a dark soliton: this solitonvortex mixed state has a dark soliton stripe nodal line plus two additional lines over each of which three vortices reside. Similar states with 8, 10 etc. vortices beside the soliton also exist, arising through subsequent bifurcations.
To corroborate the above theoretical predictions we show in Figs. 6 and 7 (again for ) some prototypical examples of states that emerge from the third excited linear branches and their corresponding BdG stability spectra. The 9vortex “crystal” is stable close to the linear limit in the sense that its BdG spectrum possesses no imaginary mode. However there is a complex mode inducing an oscillatory instability. At the vortex plus ring soliton state bifurcates from the 9vortex “crystal” inducing a complex mode. However, the latter is small and vanishes again at . The BdG spectrum of the vortex plus ring soliton state contains no purely imaginary mode but many complex modes. Note that the eigenvalue spectrum looks fairly similar to the ring dark soliton case, but the imaginary part of the modes is not due to imaginary modes which are created by bifurcations. In this case, modes with positive energy cross zero and thus obtain a negative energy. These negative energy modes then collide with modes with positive energy and create the complex modes. Moreover, we show the threesolitonstripe, as well as the mixed solitonsixvortex cluster state emerging from its first bifurcation event in Fig. 6. The threesoliton state is stable near the linear limit (of ) but becomes destabilized at (and then further so at , and ). The first bifurcation gives rise to the very weakly unstable solitonvortex state predicted above; the critical point for this bifurcation is found to be at , once again in very good agreement with the full numerical result. The present approach can naturally be extended to a multitude of additional states which, however, are expected to be dynamically unstable; therefore, having presented the most fundamental ones, we will not proceed further with such considerations.
4 Conclusions and Future Directions
In conclusion, in the present work we have shown that a detailed understanding of emergent vortex cluster states is possible through a nearlinearlimit approach. This involves identifying the possible linear states and tracking systematically the symmetrybreaking bifurcations that can arise from them. This allowed us not only to discuss a cascade of bifurcations from the first excited state (in the form of aligned vortex clusters), but to also reveal a broad class of states emerging from the second and thirdexcited states. These were not only rectilinear states (with “solitontype” stripes), but also soliton rings and rings with vortices, vortex quadrupoles or ninevortexcrystals, as well as various clusters of vortices that derive from some of these states, including vortex hexagons, octagons, decagons, states (of vortices sitting along stripes), solitonvortex states, and so on.
We were also able, based on the general bifurcation structure of the problem, to reveal which ones among these states are expected to be most robust, such as the vortex dipole or quadrupole, and some of the emerging vortex clusters, such as the or the solitonvortex state.
It would be especially interesting to extend this picture to different dimensions. On the one hand, one can consider effects of anisotropy (by changing the strengths of the two different trapping directions). This should enable a “dimensiontranscending” picture, as extreme anisotropies should allow to observe how the system transitions between quasione and genuinely twodimensional dynamics. On the other hand, it would be particularly relevant and interesting to attempt to extend such considerations into threedimensional settings, and understand how relevant ideas generalize and potentially give rise to structures such as vortex rings that have been observed in pertinent experiments [33, 42]. Such investigations are currently in progress.
Acknowledgments
P.G.K. gratefully acknowledges support from NSFDMS0349023 (CAREER), NSFDMS0806762 and the Alexander von Humboldt Foundation. The work of D.J.F. was partially supported by the Special Account for Research Grants of the University of Athens. R.C.G. acknowledges support from NSFDMS0806762. The authors also gratefully acknowledge David Hall for numerous enlightening discussions on experimental aspects of this theme.
References
 [1] L. M. Pismen, Vortices in Nonlinear Fields (Oxford Science Publications, Oxford, 1999).
 [2] R. J. Donnelly, Quantized Vortices in Helium II (Cambridge University Press, New York, 1991); D. R. Tilley and J. Tilley, Superfluidity and Superconductivity (IOP Publishing, Philadelphia, 1990).
 [3] A. J. Chorin and J. E. Marsden, A Mathematical Introduction to Fluid Mechanics (SpringerVerlag, New York, 1993).
 [4] P. G. Kevrekidis, D. J. Frantzeskakis, and R. CarreteroGonzález, Emergent Nonlinear Phenomena in BoseEinstein Condensates (SpringerVerlag, Berlin, 2008).
 [5] A. S. Desyatnikov, Yu. S. Kivshar, and L. Torner, Progress in Optics 47, 291 (2005).
 [6] M. R. Matthews, B. P. Anderson, P. C. Haljan, D. S. Hall, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. 83, 2498 (1999).
 [7] J. E. Williams and M. J. Holland, Nature 401, 568 (1999).
 [8] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard, Phys. Rev. Lett. 84, 806 (2000).
 [9] A. Recati, F. Zambelli, and S. Stringari, Phys. Rev. Lett. 86, 377 (2001).
 [10] S. Sinha and Y. Castin, Phys. Rev. Lett. 87, 190402 (2001).
 [11] I. Corro, R. G. Scott, and A. M. Martin, Phys. Rev. A 80, 033609 (2009).
 [12] K. W. Madison, F. Chevy, V. Bretin, and J. Dalibard, Phys. Rev. Lett. 86, 4443 (2001).
 [13] C. Raman, J. R. AboShaeer, J. M. Vogels, K. Xu, and W. Ketterle, Phys. Rev. Lett. 87, 210402 (2001).
 [14] R. Onofrio, C. Raman, J. M. Vogels, J. R. AboShaeer, A. P. Chikkatur, and W. Ketterle, Phys. Rev. Lett. 85, 2228 (2000).
 [15] D. R. Scherer, C. N. Weiler, T. W. Neely, and B. P. Anderson, Phys. Rev. Lett. 98, 110402 (2007).
 [16] A.E. Leanhardt, A. Görlitz, A. P. Chikkatur, D. Kielpinski, Y. Shin, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett. 89, 190403 (2002); Y. Shin, M. Saba, M. Vengalattore, T. A. Pasquini, C. Sanner, A. E. Leanhardt, M. Prentiss, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett. 93, 160406 (2004).
 [17] T. W. Neely, E. C. Samson, A. S. Bradley, M. J. Davis, and B. P. Anderson, Phys. Rev. Lett. 104, 160401 (2010).
 [18] D. V. Freilich, D. M. Bianchi, A. M. Kaufman, T. K. Langin, and D. S. Hall, Science 329, 1182 (2010).
 [19] J. A. Seman, E. A. L. Henn, M. Haque, R. F. Shiozaki, E. R. F. Ramos, M. Caracanhas, P. Castilho, C. Castelo Branco, K. M. F. Magalhes, and V. S. Bagnato, Phys. Rev. A 82, 033616 (2010).
 [20] T. W. B. Kibble, J. Phys. A 9, 1387 (1976); W. H. Zurek, Nature 317, 505 (1985); W.H. Zurek, Phys. Rep. 276, 177 (1996).
 [21] C. N. Weiler, T. W. Neely, D. R. Scherer, A. S. Bradley, M. J. Davis, B. P. Anderson, Nature 455, 948 (2008).
 [22] P. Kuopanportti, J.A.M. Huhtamäki and M. Möttönen, arXiv:1011.1661.
 [23] L.C. Crasovan, V. Vekslerchik, V. M. PérezGarcía, J. P. Torres, D. Mihalache, and L. Torner, Phys. Rev. A 68, 063609 (2003).
 [24] M. Möttönen, S. M. M. Virtanen, T. Isoshima, and M. M. Salomaa, Phys. Rev. A 71, 033626 (2005).
 [25] V. Pietilä, M. Möttönen, T. Isoshima, J. A. M. Huhtamäki and S. M. M. Virtanen, Phys. Rev. A 74, 023603 (2006).
 [26] A. Klein, D. Jaksch, Y. Zhang, and W. Bao, Phys. Rev. A 76, 043602 (2007).
 [27] W. Li, M. Haque and S. Komineas, Phys. Rev. A 77, 053610 (2008).
 [28] J. Brand and W. P. Reinhardt, Phys. Rev. A 65, 043612 (2002).
 [29] R. S. MacKay, in Hamiltonian Dynamical Systems, edited by R. S. MacKay and J. Meiss (Hilger, Bristol, 1987), p.137.
 [30] D. J. Frantzeskakis, J. Phys. A: Math. Theor. 43, 213001 (2010).
 [31] E. A. Kuznetsov and S. K. Turitsyn, Zh. Eksp. Teor. Fiz. 94, 119 (1988) [Sov. Phys. JEPT 67, 1583 (1988)].
 [32] D. L. Feder, M. S. Pindzola, L. A. Collins, B. I. Schneider, and C. W. Clark, Phys. Rev. A 62, 053606 (2000).
 [33] 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).
 [34] S. Middelkamp, P. G. Kevrekidis, D. J. Frantzeskakis, R. CarreteroGonzález, and P. Schmelcher, Phys. Rev. A 82, 013646 (2010).
 [35] G. Theocharis, P. G. Kevrekidis, D. J. Frantzeskakis, and P. Schmelcher, Phys. Rev. E 74, 056608 (2006).
 [36] T. Zibold, E. Nicklas, C. Gross and M.K. Oberthaler, Phys. Rev. Lett. 105, 204101 (2010).
 [37] T. Kapitula, P. G. Kevrekidis and R. CarreteroGonzález, Physica D 233, 112 (2007).
 [38] G. Theocharis, D. J. Frantzeskakis, P. G. Kevrekidis, B. A. Malomed, and Yu. S. Kivshar Phys. Rev. Lett. 90, 120403 (2003).
 [39] G. Theocharis, P. Schmelcher, M. K. Oberthaler, P. G. Kevrekidis, and D. J. Frantzeskakis, Phys. Rev. A 72, 023609 (2005); L. D. Carr and C. W. Clark, Phys. Rev. A 74, 043613 (2006).
 [40] G. Herring, L. D. Carr, R. CarreteroGonzález, P. G. Kevrekidis, and D.J. Frantzeskakis, Phys. Rev. A 77, 023625 (2008).
 [41] P.G. Kevrekidis and D.E. Pelinovsky, Phys. Rev. A 81, 023627 (2010).
 [42] N. S. Ginsberg, J. Brand and L. V. Hau, Phys. Rev. Lett. 94, 040403 (2005).