# A two-dimensional paradigm for symmetry breaking:

the nonlinear Schrödinger equation with a four-well potential

###### Abstract

We study the existence and stability of localized states in the two-dimensional (2D) nonlinear Schrödinger (NLS)/Gross-Pitaevskii equation with a symmetric four-well potential. Using a four-mode approximation, we are able to trace the parametric evolution of the trapped stationary modes, starting from the corresponding linear limits, and thus derive the complete bifurcation diagram for the families of these stationary modes. The predictions based on the four-mode decomposition are found to be in good agreement with the numerical results obtained from the NLS equation. Actually, the stability properties coincide with those suggested by the corresponding discrete model in the large-amplitude limit. The dynamics of the unstable modes is explored by means of direct simulations. Finally, while we present the full analysis for the case of the focusing nonlinearity, the bifurcation diagram for the defocusing case is briefly considered too.

## I Introduction

In the recent years, there has been a considerable effort on experimental and theoretical studies of Bose-Einstein condensates (BECs) book1 (); book2 (). Many of these studies were focused on macroscopic nonlinear structures that arise in BECs, which often have counterparts in nonlinear optics agra (). One of the appealing features of this setting is the possibility to tailor the desirable geometry of magnetic, optical, or combined traps that confine the ultracold bosonic atoms. For this reason, the analysis of the existence, stability and dynamical properties of nonlinear localized states trapped in these geometries have become a focal point of research. The theoretical analysis is enabled by the fact that a very accurate description of dilute atomic BECs is furnished, in the mean-field approximation, by the Gross-Pitaevskii (GP) equation, which is a variant of the nonlinear Schrödinger (NLS) equation. The cubic nonlinearity in the GP equation originates from the interatomic interactions, accounted for through an effective mean-field. The NLS equation in this, as well as in somewhat different forms, is relevant to a variety of alternative physical applications in nonlinear optics and other areas agra (); ourbook (); sulem ().

Among the trapping configurations available in current BEC experiments, one that has drawn considerable attention is the double-well potential (DWP). Its prototypical realization is provided by the a strong parabolic (harmonic) trap combined with a periodic potential, which can be created as an “optical lattice”, by a set of coherent laser beams illuminating the condensate book1 (); book2 (); ourbook (). The DWP was realized experimentally in markus1 (), using the magnetic field to induce the parabolic trap. The experiments reported in Ref. markus1 () revealed a variety of fundamental effects, including tunneling and Josephson oscillations for a small number of atoms, or macroscopic quantum self-trapping leading to a stable asymmetric partition of atoms between the wells for a sufficiently large number of atoms. DWPs have also inspired theoretical studies of various topics, such as finite-mode reductions, finding exact analytical results for specially designed shapes of the potential, quantum effects smerzi (); kiv2 (); mahmud (); bam (); Bergeman_2mode (); infeld (); todd (); theo (); carr (), and a nonlinear DWP (alias double-well pseudopotential), which is induced by the respective spatial modulation of the nonlinearity coefficient pseudo (). It is relevant to mention that DWPs have also been studied in the context of nonlinear optics, including twin-core self-guided laser beams in Kerr media HaeltermannPRL02 () and optically-induced waveguiding structures in photorefractive crystals zhigang ().

The aim of the present work is to extend the analysis of the DWP to a two-dimensional (2D) setting. Unlike previous studies of the trapping of quasi-2D BECs under the combined action of harmonic traps and optical-lattice potentials kody (), we implement a Galerkin-type few-mode reduction to deduce a discrete model, based on a symmetric set of four wells, which is the most natural configuration in the 2D case. The same model can also be realized in optics, using a bulk nonlinear medium with a set of four embedded waveguiding rods. In that setting, we use numerical methods to generate a bifurcation diagram for possible stationary states of the system. It is worth noting that all the states that are expected on the basis of a four-site discrete nonlinear Schrödinger (DNLS) reduction wannier () are also obtained in the continuum model with the combined parabolic and periodic potential considered herein. Furthermore, their stability, in the large nonlinearity limit, coincides with what is expected from the DNLS model.

The paper is structured as follows. In section II, we present the model and the derivation of the four-mode approximation. Numerical results are reported in section III. We present complete bifurcation diagrams of the possible stationary states for both the underlying GP equation and for its four-mode reduction. Comparison between them demonstrates very good agreement. In addition to the study of the existence and stability, evolution of unstable modes is explored too by means of direct numerical simulations. Finally, we summarize our findings in section IV, where we also discuss possible directions for further work.

## Ii The model and the Galerkin approximation

We start by presenting the basic model in the quasi-2D setting, namely the NLS/GP equation in -dimensions, which is expressed in the following dimensionless form book1 (); book2 (); ourbook ():

(1) |

where is the normalized wave function, the chemical potential ( is the propagation constant in the optical realization of the model), corresponds, respectively, to repulsive or attractive interatomic interactions in BEC (alias self-defocusing or self-focusing Kerr nonlinearity, in terms of nonlinear optics), and is the single-particle operator given by:

(2) |

In Eq. (2), is the Laplacian, while is the four-well potential, assumed to take the following form,

(3) |

with . It is clear that is composed of a harmonic trap of strength and a periodic (OL) potential with strength and period . In the following analysis below, we adopt the following representative values for parameters of the potential: , and , in which case the four smallest eigenvalues of operator are found to be

(4) |

In a weakly nonlinear setting, we implement the natural possibility of a four-mode approximation, based on a Galerkin-type expansion of and truncation of the higher-order modes. We denote the ground and first three excited eigenstates of the linear operator , shown in Fig. 1, as and . This set constitutes a natural minimum basis for the Galerkin approximation in the system of four potential wells coupled by tunneling. Eigenstates () can be chosen to be real, given the Hermitian nature of the operator . Then, solutions of Eq. (1) for values of the chemical potential in a vicinity of linear eigenvalues (4), may be approximated by linear combinations of the four eigenfuctions.

Actually, it is more convenient to use a transformed basis, {, , , }, as shown in Fig. 2, which is based on populations of the four wells. This basis is generated from the original set by a linear transformation:

(5) |

where the appropriate transformation matrix is

(6) |

Each mode () is localized in one of the four wells, with the four of them constituting an orthonormal set.

Using the new basis, we can readily reformulate the four-mode decomposition as

(7) |

with time-dependent complex amplitude , . Substituting Eq. (7) into Eq. (1) and projecting onto the orthonormal basis {, , , }, we derive, by means of straightforward algebra, the following system of four ordinary differential equations (ODEs),

(8) |

with the summation performed over . To cast these equations in a more compact form, we have defined

(9) |

where , , , and . Notice that, for the underlying eigenvalues (4), one has and . Furthermore, Eqs. (8) involve nonlinear coefficients given by overlap integrals, viz., , , , , , with ; these indices must be mutually different wherever they appear in the coefficients before the nonlinear terms.

For our choice of the parameters of the potential, the overlapping between modes is weak (see Fig. 2), therefore all the overlap integrals are much smaller than the ’s. Neglecting these small overlap terms leads to the following simplification of Eq. (8):

(10) |

It has been checked that the latter reduction of the four-mode equations very slightly affects the accuracy of the solutions, while it renders the identification of various bifurcation branches significantly easier. Furthermore, this reduction is more convenient in simulations as we may use these solutions as inputs for generating numerical solutions of the full GP system, as explained below.

In this 2D setting, we seek both real and complex stationary solutions to the ODE system. Substituting into Eq. (10), we split them into real equations for and :

(11) | ||||

(12) | ||||

(13) | ||||

(14) |

with the equations for and obtained by interchanging the indices, and , except for in and .

Looking for solutions with constant and which are integer multiples of , we reduce Eqs. (11) - (14) to a set of four algebraic equations for , which can be used to derive a complete set of stationary solutions of the four-mode truncation. These were further used as initial guesses to find numerical solutions of the full system of the GP equations. Moreover, our analysis of the four-mode system indicates that nontrivial complex solutions in this setting are only possible in the form of discrete vortices, i.e., solutions with phase sets , peli2d (), which have been studied in detail in Refs. kody () and todd_chaos () (we also briefly consider them here).

## Iii Numerical results

### iii.1 Attractive interactions

Let us first present results of numerical simulations pertaining to attractive interactions (alias self-focusing nonlinearity) case, i.e., in Eq.(1). Our basic bifurcation diagram, shown in Fig. 3, displays the squared norm of the solution (which physically describes the number of atoms in BECs or the power in optics), , as a function of the chemical potential . The left panel of Fig. 3 presents the full numerical bifurcation diagram, which involves twelve real and one complex solutions (for the latter solution, is the same as one of the real branches, hence this branch is not visible as a separate curve in the diagram). The companion diagram in the right panel is obtained from the above-mentioned algebraic system for stationary solutions produced by the the four-mode reduction, and demonstrates good agreement with its numerical counterpart.

The twelve real branches are labeled mainly according to their relation to the populations of the four wells. To support our explanation, we introduce a symbolic representation that we developed in the form of matrices, labeling different waveforms that arise in the diagram, as follows: , , , , , , , , , , , and . In this representation, , and have the obvious meaning, by indicating that a particular well is or is not populated, and the phase of the wave function in it being or in the cases of and , respectively, when populated. Symbol denotes either a small (but nonzero) population in one of the wells, or a symmetry-breaking effect (when some of the density peaks feature values , thus being slightly different from ). The labeling is then defined as follows: branches A1-A4 have the same amplitude at the four wells as long as they are populated, branches B1-B3 feature two pairs of peaks with different amplitudes, branches C1-C4 have three different amplitudes, while D1 has all of its four peaks different. The waveforms in the top rows of Figs. 4-7 display prototypical realizations of the relevant branches. Their stability properties are illustrated, as a function of eigenvalue parameter , in the bottom rows.

We will now explain in detail solutions appearing in the full bifurcation diagram, starting from the linear limits (). First, we look at the group of solutions related to branch A1, as shown in the bottom left panel of Fig. 3. This branch arises from the symmetric linear mode at , i.e., the ground state in the linear limit, . Accordingly, A1 features four identically populated wells. The analysis demonstrates that it is stable near the linear limit, but is soon destabilized, due to the emergence of branches C1 and B1, through subcritical and supercritical pitchfork bifurcations, respectively, around . In other words, there are two consecutive steady-state bifurcations, in the language of Ref. dellnitz (), in two different subspaces, in which one unstable solution, C1, collides with A1 and, simultaneously, a pair of eigenvalues emerges on the real axis for the A1 branch with the decrease of (in a subcritical pitchfork); then, a super-critical pitchfork takes place in another subspace, in which B1 retains only one real pair, while another pair passes through the origin along the A1 branch. The actual “pitchfork” cannot be visualized here in the usual manner, because any of the four equivalent versions of B1 (obtained by rotation) have the same value of , and so are represented by the same branch in the graph. Branch B1, which is unstable due to a pair of real eigenvalues in the linearization around it throughout its domain of existence, features two of the wells on one side being less populated than the other two. Configuration B1 becomes increasingly more asymmetric as it deviates from A1. A noteworthy feature is shown by branch C1, which bifurcates from A1 at almost the same place as B1: after having emerged, it tends to be located on the left of A1 as are all other branches bifurcating from A1. However, within a narrow interval of , its norm decreases () slightly before starting to grow as usual (). Naturally, when the norm decreases the solution is destabilized and then it remains stable after the turning point, which is explained by the well-known Vakhitov-Kolokolov criterion vk (). Branch A1 is endowed with two identical pairs of real eigenvalues by B1 and C1 upon their bifurcation (which is shown as the dashed line in the bottom left panel of Fig. 4.) As grows further, a subsequent bifurcation, at , leading to the emergence of branch B2, adds yet another real eigenvalue pair to A1; this means that A1 possesses in total three real eigenvalue pairs for sufficiently large vaues of . Branch B2 features two wells on the diagonal which are less populated than the other two, and it is unstable, with two pairs of real eigenvalues, near the point where it is generated by the bifurcation from A1; however, one of the pairs is eliminated by the emergence of a new branch, C2, from B2 through a pitchfork shortly afterwards. Branch B2 then remains unstable with one real pair, while C2 (with three principal sites, one of which is of lesser amplitude than the other two) is unstable due to two real eigenvalue pairs.

This description encompasses all branches of stationary solutions which can be traced back to the ground state of the linear system. Detailed information for the wave function profiles and the development of the real eigenvalues associated to them is presented in Figs. 4-6.

Next we turn to the states originating from the second linear mode, as shown in the bottom right panel of Fig. 3. Branch A2 starts from the respective eigenvalue, , which pertains to the first and second (degenerate at the linear limit) excited states. This branch emerges as an unstable one, carrying a real eigenvalue pair. The respective wave function profile features four wells populated with the same amplitude but out-of-phase between the two sides, see Fig. 4. Branch B3 emerges from A2 through a supercritical pitchfork at , lending another real eigenvalue pair to A2. Similar to the case of B1 (as it separates from A1), in the B3 state two wells on the one side tend to be less populated than the other two, as this branch moves further from A2. Branch B3 remains unstable with one real eigenvalue pair, until getting stabilized by another pitchfork bifurcation which takes place at ; this simultaneously gives rise to a new unstable branch, D1, that has different populations in all four wells. Notice that both B3 and D1 pass through a Hamiltonian-Hopf bifurcation (alias 1:1 resonance, in terms of Ref. dellnitz ()), which means that, in the relevant parametric interval ( for B3; for D1), B3 is destabilized by a complex quartet of small eigenvalues in the linearization around the stationary solution, while D1 remains unstable, but with one real eigenvalue pair and a complex quartet (the dashed-dotted lines in the bottom right panels of Figs. 5 and 7 refer to this effect).

Furthermore, branch A4 bifurcates from the same linear mode as A2. It is unstable near the linear limit due to a Hamiltonian-Hopf bifurcation, but with the increase of it becomes stable. Branch A4 features a waveform in which only two wells lying on the diagonal are populated, with the same amplitude but opposite signs.

Branch A3 arises from the third excited linear mode at . In this case, two wells on the diagonal are populated with equal amplitudes, while in the other two the amplitudes are of opposite signs. It is the unique stationary solution which remains stable across the entire bifurcation diagram (despite the fact that it has three pairs of purely imaginary eigenvalues with negative Krein signature sandstede (), which in principle, can give rise to Hamiltonian-Hopf bifurcations).

Finally, branches C3 and C4, which are located slightly below A2 in Fig. 3, correspond to a pair of states which arise through a saddle-node bifurcation at some critical value of chemical potential (, for the model’s parameters chosen in the present case.) Branch C4 (the one with higher values of ) is unstable with a real eigenvalue pair, while C3 remains stable, except inside a short instability interval, which is accounted for by a Hamiltonian-Hopf bifurcation.

In addition to the above real stationary states, we have also found complex solutions in the form of vortices kody (); todd_chaos (). A typical example of such a solution is shown in Fig. 8. Throughout the regime of parameters considered herein, such solutions have been found to be linearly stable.

It is interesting to note that, for all the solutions considered herein, in the large- limit their stability characteristics coincide with what can be suggested by the DNLS model considered in Ref. peli2d () (see also Refs. peli1d (); peli3d () for the corresponding 1D and 3D stability results). Gross features of these findings are that, whenever two adjacent sites are in-phase, a real eigenvalue pair is expected to emerge due to their interaction, while whenever such sites are out-of-phase, the relevant eigenvalue is expected to be imaginary multip (), but with negative Krein signature peli2d (), which implies a potential for a Hamiltonian-Hopf bifurcation. It should also be noted that, in the limit of the infinite lattice, it is naturally expected that the asymmetries observed herein in many of the branches will disappear (i.e., the amplitudes in different wells will be equal) – see also a relevant discussion in Ref. 3well ().

### iii.2 Repulsive interactions

We now briefly discuss the case of repulsive interactions (alias self-defocusing nonlinearity), corresponding to in Eq. (1), with an objective to highlight its similarities with and differences from the case of attractive interactions. The bifurcation diagram for the model is displayed in Fig. 9.

The solutions are labeled so as to match the self-focusing case, by means of the appropriate staggering transformation hadistuff (). The latter effectively converts the defocusing nonlinearity into a focusing one by changing the relative phase of nearest-neighbors from to and vice versa, while preserving the relative phase of next-nearest-neighbors. In this way, each solution in the defocusing case is linked to its counterpart in the focusing model through this transformation. Following this relation, and adopting the same matrix symbolic representation used for the focusing case in section III.1, the branches of solutions are labeled as follows: , , , , , , , , , , , and .

Thus, in this case, the symmetric ground state of the system is A3, which is stable for arbitrary values of . Branch A2 is immediately unstable, starting from the linear limit. B3 bifurcates from A2 and remains unstable before getting stabilized through giving birth to D1 (and then becoming destabilized again). Branch A1 is stable near the linear limit, but is subsequently destabilized due to bifurcations that give rise to B1 and C1, and an additional real eigenvalue pair arises at higher value of due to the emergence of B2, from which another new branch, namely C2, arises in turn. Branches C3 and C4 exist for a while (when is large enough), and then collide at (for the values of parameters adopted herein). The types of the bifurcations, the emergence of the corresponding solutions, and the corresponding stability properties were found to be in direct correspondence to the case of self-focusing nonlinearity, provided that one takes into account the staggering transformation relating the repulsive and attractive case models as indicated above.

### iii.3 Dynamics

We now proceed to investigate the evolution of unstable states in the model with the self-focusing nonlinearity. To this end, for each unstable branch, a small perturbation is added to the most unstable eigendirection of the linearization near the original stationary solution, at . Results of the simulations are presented in Fig. 10.

Panel (a) shows the behavior of solution A1, which, as a result of the instability, starts oscillating between a state where all four wells are populated and one in which only two diagonal wells are not empty. Panel (b) depicts a periodic oscillatory behavior of A2, also between four and two populated sites, but in this case the continuously populated sites are adjacent to each other. Unstable mode A4 (panel (c)) features only two non-empty wells, with the symmetry-breaking instability resulting in the enhanced population of one of the two. As explained in section III.1, modes B1, B2 and B3 have two very weakly populated wells, in comparison with the other two. Since we employ isosurface (where is the half of the maximum density at ) to plot the dynamics in Fig. 10, the evolution in the weakly populated wells is not visible (which indicates that they play a minor role in the dynamics). Panel (d) shows that mode B1 sustains a symmetry breaking similar to that of A4, but between adjacent sites. On the other hand, modes B2 and B3 appear to be oscillating between the two dominant wells roughly periodically, as shown in panels (e) and (f). Mode C2 [in panel (g)] oscillates between three and two populated sites (the seemingly empty well is actually a weakly populated one, similarly to the modes of type B, see above). Mode C3, whose weak instability is caused by a quartet of eigenvalues, is also “breathing” within the respective set of three predominantly populated wells, as shown in panel (h). Finally, mode C4 [shown in panel (i)] involves a complex symmetry-breaking pattern, with different numbers of wells populated at different times, while mode D1 (panel (j)) oscillates between three- and two-well asymmetric configurations.

## Iv Conclusion

In this work, we have studied stationary and dynamical properties of the two-dimensional nonlinear Schrödinger/Gross-Pitaevskii equation, including a four-well external potential, with both signs of the nonlinearity, self-attractive (focusing) and self-repulsive (defocusing). The model applies to a BEC confined in a highly anisotropic harmonic trap, with the transverse confining frequency being much smaller than the axial one, which results in a planar configuration. In this context, the four-well potential can be generated by a combination of the transverse part of the harmonic trap and an optical lattice. The same model may describe the propagation of an optical beam in a bulk nonlinear medium with an embedded four-channel guiding structure.

In our analysis, first we developed a four-mode approximation, which strongly simplifies the identification of stationary solutions. Using this approximation, we were able to find the four (two of which are identical) symmetric and antisymmetric linear modes, and all branches of asymmetric solutions emerging from them in the model with the self-focusing nonlinearity. The linear-stability analysis demonstrated how pitchfork and saddle-node bifurcations change the stability of the branches. We have shown that in the limit of strong nonlinearity, properties of localized modes in the model with either sign of the nonlinearity can be, roughly, understood on the basis of earlier results pertaining to the corresponding discrete NLS model. We have also described the evolution of all unstable solutions, observing, typically, the emergence of symmetry-breaking instabilities and the emergence of respective oscillating solutions.

It would be interesting to investigate how these four-site configurations may be embedded into a larger potential pattern, with or wells, and examine whether the symmetry-breaking bifurcations considered above are sustained (or how they are modified) within the larger pattern. In this context, a conjecture that would be worthwhile proving is that, in an infinite periodic lattice formed by potential wells, the nonlinearity can support 2D solitons and localized vortices with various symmetries, but not confined asymmetric states. This conjecture is suggested by results reported for infinite linear infinite () and nonlinear Barcelona () potential lattices.

The work of D.J.F. was partially supported by the Special Research Account of the University of Athens. P.G.K. gratefully acknowledges support from NSF-DMS-0349023 and NSF-DMS-0806762, as well as from the Alexander von Humboldt Foundation.

## References

- (1) L. P. Pitaevskii, 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) Yu. S. Kivshar and G. P. Agrawal, Optical solitons: from fibers to photonic crystals (Academic Press, San Diego, 2003).
- (4) P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero-González (eds.), Emergent nonlinear phenomena in Bose-Einstein condensates. Theory and experiment (Springer-Verlag, Berlin, 2008).
- (5) C. Sulem and P. L. Sulem, The Nonlinear Schrödinger Equation (Springer-Verlag, New York, 1999).
- (6) M. Albiez, R. Gati, J. Fölling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, Phys. Rev. Lett. 95, 010402 (2005).
- (7) S. Raghavan, A. Smerzi, S. Fantoni, and S. R. Shenoy, Phys. Rev. A 59, 620 (1999); S. Raghavan, A. Smerzi, and V. M. Kenkre, Phys. Rev. A 60, R1787 (1999); A. Smerzi and S. Raghavan, Phys. Rev. A 61, 063601 (2000).
- (8) E. A. Ostrovskaya, Yu. S. Kivshar, M. Lisak, B. Hall, F. Cattani, and D. Anderson, Phys. Rev. A 61, 031601(R) (2000).
- (9) K. W. Mahmud, J. N. Kutz, and W. P. Reinhardt, Phys. Rev. A 66, 063607 (2002).
- (10) V. S. Shchesnovich, B. A. Malomed, and R. A. Kraenkel, Physica D 188, 213 (2004).
- (11) D. Ananikian and T. Bergeman, Phys. Rev. A 73, 013604 (2006).
- (12) P. Ziń, E. Infeld, M. Matuszewski, G. Rowlands, and M. Trippenbach, Phys. Rev. A 73, 022105 (2006).
- (13) T. Kapitula and P. G. Kevrekidis, Nonlinearity 18, 2491 (2005).
- (14) G. Theocharis, P. G. Kevrekidis, D. J. Frantzeskakis, and P. Schmelcher, Phys. Rev. E 74, 056608 (2006).
- (15) D. R. Dounas-Frazer, A. M. Hermundstad, and L. D. Carr, Phys. Rev. Lett. 99, 200402 (2007).
- (16) T. Mayteevarunyoo, B. A. Malomed, and G. Dong. Phys. Rev. A 78, 053601 (2008).
- (17) C. Cambournac, T. Sylvestre, H. Maillotte , B. Vanderlinden, P. Kockaert, Ph. Emplit, and M. Haelterman, Phys. Rev. Lett. 89, 083901 (2002).
- (18) P. G. Kevrekidis, Z. Chen, B. A. Malomed, D. J. Frantzeskakis, and M. I. Weinstein, Phys. Lett. A 340, 275 (2005).
- (19) K. J. H. Law, L. Qiao, P. G. Kevrekidis, and I. G. Kevrekidis, Phys. Rev. A 77, 053612 (2008); K. J. H. Law, P. G. Kevrekidis, B. P. Anderson, R. Carretero-González, and D. J. Frantzeskakis, J. Phys. B: At. Mol. Opt. Phys. 41, 195303 (2008).
- (20) G. L. Alfimov, P. G. Kevrekidis, V. V. Konotop, and M. Salerno, Phys. Rev. E 66, 046608 (2002).
- (21) T. Kapitula, P. G. Kevrekidis, and D. J. Frantzeskakis, Chaos 18, 023101 (2008).
- (22) D. E. Pelinovsky, P. G. Kevrekidis, and D. J. Frantzeskakis, Physica D 212, 20 (2005).
- (23) N. G. Vakhitov and A. A. Kolokolov, Radiophys. Quantum Electron. 16, 783 (1973).
- (24) D. E. Pelinovsky, P. G. Kevrekidis, and D. J. Frantzeskakis, Physica D 212, 1 (2005).
- (25) M. Lukas, D. E. Pelinovsky, and P. G. Kevrekidis, Physica D 237, 339 (2008).
- (26) T. Kapitula, P. G. Kevrekidis, and B. A. Malomed, Phys. Rev. E 63, 036604 (2001).
- (27) T. Kapitula, P. G. Kevrekidis, and Z. Chen, SIAM J. Appl. Dyn. Sys. 5, 598 (2006).
- (28) P. G. Kevrekidis, H. Susanto, and Z. Chen, Phys. Rev. E 74, 066606 (2006).
- (29) M. Dellnitz, I. Melbourne, and J. E. Marsden, Nonlinearity 5, 979 (1992).
- (30) T. Kapitula, P.G. Kevrekidis and B. Sandstede, Physica D 195, 263 (2004).
- (31) R. Driben, B. A. Malomed, A. Gubeskys, and J. Zyss, Phys. Rev. E 76, 066604 (2007); R. Driben and B. A. Malomed, Eur. Phys. J. D 50, 317 (2008).
- (32) Y. V. Kartashov, B. A. Malomed, V. A. Vysloukh, and L. Torner, Opt. Lett. 34, 770 (2009).