# Surface solitons in trilete lattices

###### Abstract

Fundamental solitons pinned to the interface between three semi-infinite one-dimensional nonlinear dynamical chains, coupled at a single site, are investigated. The light propagation in the respective system with the self-attractive on-site cubic nonlinearity, which can be implemented as an array of nonlinear optical waveguides, is modeled by the system of three discrete nonlinear Schrödinger equations. The formation, stability and dynamics of symmetric and asymmetric fundamental solitons centered at the interface are investigated analytically by means of the variational approximation (VA) and in a numerical form. The VA predicts that two asymmetric and two antisymmetric branches exist in the entire parameter space, while four asymmetric modes and the symmetric one can be found below some critical value of the inter-lattice coupling parameter – actually, past the symmetry-breaking bifurcation. At this bifurcation point, the symmetric branch is destabilized and two new asymmetric soliton branches appear, one stable and the other unstable. In this area, the antisymmetric branch changes its character, getting stabilized against oscillatory perturbations. In direct simulations, unstable symmetric modes radiate a part of their power, staying trapped around the interface. Highly unstable asymmetric modes transform into localized breathers traveling from the interface region across the lattice without significant power loss.

###### keywords:

soliton complex, spontaneous symmetry breaking^{†}

^{†}journal: Physica D

## 1 Introduction

Surface modes, which are a special type of waves localized at interfaces between different media, were first predicted as localized Tamm electronic states at the edge of a truncated periodic potential rad1 (). In optics, it was predicted theoretically and confirmed experimentally that the nonlinear self-trapping of light near the edge of a waveguide array with the self-focusing nonlinearity can lead to the formation of discrete surface solitons rad2 (); rad3 (). Various settings for the creation of surface solitons were also proposed for Bose-Einstein condensates (BECs) bec (). The general framework for the description of such localized patterns is provided by the discrete nonlinear-Schrödinger (DNLS) equations, with appropriate boundary conditions book ().

The solitary surface modes exist above a certain threshold power and, in a certain domain of the parameter space, different surface modes can exist simultaneously. The surface modes may be understood as discrete optical solitons localized near the surface rad4 ()-rad5 (), or as lattice solitons pinned by defects rad5a ()-rady (). Surface solitons supported by truncated superlattices were investigated too He-super (). Moreover, the light localization in self-defocusing nonlinear media in the form of surface gap solitons has been predicted and observed in Refs. rad7 (); rad8 (); rad8a (); rad6 (), and the concept of multi-gap surface solitons, i.e., mutually trapped surface modes with components associated with different spectral gaps, was put forward rad6 (); rad6a () (multi-gap, alias ”inter-gap”, or ”semi-gap”, solitons are also known in uniform lattice media inter ()). A short review of surface solitons in discrete systems was given in Ref. He ().

The studies of surface modes have shown that nonlinear discrete photonic and matter-wave systems support spatially localized states with sundry symmetries (which can be controlled by the insertion of suitable defects into the system) rad5d (); bec (). Related to this is the possibility of the spontaneous symmetry breaking (SSB) in symmetric dual-core systems, with a linear coupling between the two parallel cores. In fact, the SSB bifurcation, which destabilizes symmetric states and gives rise to asymmetric ones, was originally predicted in terms of the self-trapping in discrete systems 13 (). In the physically important model of dual-core nonlinear optical fibers, the SSB instability was discovered in Ref. 14 (), and the respective bifurcations for continuous-wave states were studied in detail in Ref. 15 (), for various types of the intra-core nonlinearities. Further, the SSB was studied for solitons (rather than continuous waves) in the model of the dual-core fiber with the cubic (Kerr) nonlinearity 18 (); 17 (), for gap solitons in the models of dual-core 19 () and tri-core 20 () fiber Bragg gratings, and for matter-wave solitons in the BEC loaded into a dual-core potential trap, that may be combined with a longitudinal optical lattice 26a (). In addition, the SSB was also analyzed in models describing optical media with quadratic 25 () and cubic-quintic 26 () nonlinearities.

The variational approximation (VA) Progress () makes it possible to study the SSB in dual-core systems in an analytical form 18 (); VA (). It is relevant to mention that the VA may allow one not only to describe fundamental localized modes, but also correctly predict their stability 26b ().

Continuing the investigation of the surface fundamental modes in coupled one-dimensional (1D) lattice system VA (); symnash (), in this work we study soliton complexes formed at the interface of three identical semi-infinite lattices which form a trilete configuration, see Fig. 1 below. Such lattices can be written in bulk silica by means of femtosecond optical pulses, which is the technique which has made the creation of other types of surface solitons possible Jena (). A major motivation for studying such trilete systems is that they provide a fundamental setting for the study of the SSB which is different from the well-studied dual-core configurations 20 (); Skryabin ().

The present system is modeled by three DNLS equations coupled at a single lattice site. In the framework of this system, we consider complexes built of three fundamental surface solitons, each carried by one of the three chains. Actually, these fundamental solitons are realizations of Tamm states in the present setting. We show that they form families of symmetric and asymmetric complexes, which exhibit the power-threshold behavior. Their stability and propagation dynamics are examined in detail, using the VA and numerical calculations. The existence regions for all families of the soliton complexes are produced by means of both approaches, the analytical results being quite close to the numerical ones. Stability properties, predicted by dint of the numerically implemented linear-stability analysis, are verified by direct numerical simulations of the soliton evolution.

The rest of the paper is structured as follows. The model is formulated in Section 2, which is followed by the consideration of the VA in Section 3. Numerical results are collected in Section 4, and the paper is concluded by Section 5.

## 2 The model

As said above, the trilete lattice configuration is formed by three identical semi-infinite sublattices linked at one site. The inter-site coupling constant is scaled to be inside the lattices, while a different constant, , accounts for the linkage between the lattices, having the same sign as , see Fig. 1. Thus, the triple lattice is modeled by the following system of three coupled DNLS equations,

(1) | |||||

where is the propagation distance (assuming that the chains represent three semi-infinite arrays of optical waveguides), is the discrete coordinate in the chain ( is the total number of sites in each lattice that was used in actual numerical calculations), is the Kronecker’s symbol, and rescaling is used to make the on-site self-attraction coefficient equal to .

Soliton solutions to Eqs. (1) are looked for in the usual form,

where , and are real discrete functions, and is the propagation constant. The corresponding stationary equations,

(2) | |||||

can be derived from the Lagrangian,

(3) |

where the intrinsic Lagrangians of the three semi-infinite chains are, respectively,

and the last term in Eq. (3) accounts for the coupling between them.

## 3 The variational approximation

The analytical approach to the study of the soliton solutions may be based on the VA (variational approximation) Progress (), which was adapted to discrete systems in several earlier works VA (); 26b (); symnash (); discreteVA (); mi (). Following this method, we adopt the following ansatz:

(4) |

with amplitudes , and treated as variational parameters. As concerns inverse width , following Ref. Kaup () we fix it through a solution of the linearized version of Eqs. (2) for “tails” of the discrete solitons at , from where it follows

(5) |

hence the solitons may exist for values of the propagation constant . Note that Eq. (5) yields , which is essential for the validity of analytical results displayed below – see Eqs. (8) and (10).

The substitution of ansatz (4) into Eq. (3) yields the corresponding effective Lagrangian,

(6) |

The Euler-Lagrange equations for amplitudes , and are derived from this Lagrangian in the form of

or, in the explicit form,

(7) | |||||

These equations allow us to predict the existence of different symmetric and asymmetric interface modes.

### 3.1 Existence regions for the interface solitons

The solution for symmetric solitons, with , is easily obtained from Eqs. (7):

(8) |

The dependence of this amplitude on for fixed , i.e., , is plotted in Fig. 2(a). It follows from Eq. (8) and is shown in Fig. 2(b) that the existence of the VA-predicted symmetric solutions at a given value of is limited to

(9) |

The variational calculations have shown that the SyS complexes with very small power are created near the . In this parameter region, when the nonlinear interaction is negligible, we can actually analyze the corresponding linear trilete lattice system. The straightforward analysis shows that only symmetric linear surface localized complexes with arbitrary amplitude can be formed in the linear trilete lattice system. These linear modes exist exactly at , which coincides with the critical value of the inter-lattice coupling constant () for existence of symmetric nonlinear modes. Therefore, we expect that the stable nonlinear SyS branch with small power emerges from the linear symmetric localized mode.

The critical value of at which the branch of asymmetric solitons bifurcates from the symmetric family is predicted by solving Eqs. (7) for the soliton’s amplitudes with infinitesimal differences between them: . Straightforward algebraic manipulations, which include the linearization with respect to and , yield the value of at the SSB bifurcation point: , see Fig. 3. Asymmetric solutions exist if the linear coupling is weaker than at the bifurcation point, i.e., at .

An analytical solution for asymmetric modes can be easily found in the particular case when one amplitude vanishes, . Then two other amplitudes are

(10) |

This solution may be naturally called antisymmetric. Note that, unlike symmetric solution (5), it exists for all values of . In fact, the antisymmetric solution is identical to its counterpart found in the two-chain system in Ref. mi (), although this does not mean that the stability properties of the antisymmetric solutions are identical in the two systems.

Another analytical asymmetric solution (which may be called an isosceles mode) has . In this case, can be eliminated,

(11) |

while has to be found as a numerical solution of the remaining equation,

(12) |

In accordance with the above result that determines the location of the SSB bifurcation, six different branches of solutions to Eq. (12) (three with positive and three with negative ) exist in the area of of parameter space , while only two branches (one with positive and one with negative ) are found at , as shown in Fig. 3 for . In this figure, only branches with and are shown. It is worthy to mention that, in the region of in the figure, all the asymmetric branches, with numbers 1, 2, and 3, coexist with the symmetric branch and the antisymmetric one with ; further, in the region of , only branch coexists with the same pair, while in the region of only branches and exist. It is possible to prove that solutions of Eq. (7) with all the three amplitudes different, , do not exist, which was also confirmed numerically.

Finally, it is relevant to mention that the asymptotic form of the stationary solutions can be easily understood for both and . Indeed, in the limit of the system splits into three uncoupled semi-infinite lattices. In this case, for fixed propagation constant , one may have either the usual single-component surface soliton rad2 (), or the zero solution. Accordingly, Fig. 3 demonstrates that the amplitudes of all modes at take either a fixed value, which actually corresponding to the single semi-infinite lattice, or vanish.

### 3.2 Stability of the fundamental solitons (the Vakhitov-Kolokolov criterion)

The stability of the discrete solitons predicted by the VA can be estimated by dint of the Vakhitov-Kolokolov (VK) criterion, according to which the necessary condition for the stability of fundamental solitons is , where is the total power (norm) of the soliton discreteVA (); isrl (). The total power of the solutions corresponding to ansatz (4) is

(13) |

where the relation between and is given by Eq. (5).

For the symmetric solution (8) with , dependence is plotted in Fig. 4(a) for fixed and . The slope of the curve is positive in the whole existence region of the symmetric solutions, which, according to the VK criterion, indicates the possibility for the existence of stable solutions.

The curve for the asymmetric solution branches is plotted in Fig. 4(b) for three different values of . Only the asymmetric branch which corresponds to amplitudes in Fig. 3 is VK-unstable in the narrow area of the existence region, at close to in Fig. 4(b). Other asymmetric solutions might be stable according to the VK criterion.

For the antisymmetric solutions with and , the power curves are presented in Fig. 4(c). These solutions also might be stable according to the VK criterion in the whole existence region for . For higher values of , a VK-unstable region is located near the above-mentioned existence border, [see Eq. (5)].

## 4 Numerical results

The predictions of the VA were tested by solving stationary equations (2), using an algorithm based on the modified Powell minimization method nashi (). The initial guess for constructing fundamental solitons centered at the interface of three linked semi-infinite chains (Fig. 1) was taken in the form of for symmetric modes or for asymmetric ones (with different ), while the amplitudes of the lattice field at other sites were set to be . Eventually, solutions different from symmetric ones were found solely with , or with , as predicted analytically. The results presented here are obtained for identical coupled chains of length .

The stability of the stationary modes was first checked through the linear-stability analysis. As a result, the eigenvalue (EV) spectrum for modes of small perturbations was found, following the procedure developed in Refs. symnash (); nashi (). The calculations were performed in parameter plane (, ). These results were verified by direct numerical simulations of the full system of equations (1). The simulations relied upon a numerical code which used the sixth-order Runge-Kutta algorithm, as in Refs. symnash (); nashi (). The simulations were initialized by taking stationary soliton profiles, to which random perturbations were added.

Typical shapes of symmetric and asymmetric solitons found in the numerical form are displayed in Fig. 5, well fitting the corresponding VA predictions. The respective dependencies of the solitons’ amplitudes, A and B, on coupling constant are displayed for all types of the solitons in Fig. 6 (a), while the corresponding dependencies are shown in Fig. 7. From these figures it can be concluded that both the amplitude and power characteristics give qualitatively the same information about the localized surface modes. The numerical results demonstrate that the symmetric and three asymmetric branches coexist in certain bounded regions of the parameter space, similar to what was predicted in Fig. 3.

In general, the comparison of the numerical and variational results demonstrates that the predictions of the VA for the existence region of the symmetric and asymmetric isosceles complexes of types and are very accurate. On the other hand, the numerically found existence regions for the antisymmetric solitons of type and antisymmetric complexes are bounded, on the contrary to the prediction of the VA. For example, the VA predicts that the antisymmetric complexes can exist for arbitrary and , while the numerical calculation shown the appearance of the upper limit with respect to at fixed , see Fig. 6 (b). The upper existence boundary for the antisymmetric solutions was found at extremely high values of the inter-chain coupling constant, .

The linear-stability analysis shows that a stability window exists for the symmetric soliton complexes in the region between [i.e., for close to , as it follows from Eq. (9)] and , see Fig. 8(a). These symmetric complexes are formed by solitons which are wider and possess smaller amplitudes than their counterparts belonging to unstable symmetric complexes.

The dynamics of the stable symmetric soliton complex is illustrated by Fig. 9(a) which shows the evolution of its component. In the rest of their existence region, the symmetric complexes are unstable, their instability being accounted for by purely real EV pairs. Under small perturbations, an unstable symmetric soliton complex sheds off a part of its power and relaxes into a trapped interface asymmetric breathing complex with a smaller power. The dynamics of the component belonging to the unstable symmetric complex is shown in Fig. 9(b).

The symmetric solitons with correspond to the usual surface solitons in uncoupled semi-infinite lattices. It is well known that such surface modes are stable in almost the whole existence region rad2 (), which is provided by the balance between the interaction of the soliton with the surface and the bulk lattice. The instability of the symmetric soliton complex in the coupled lattice system is related to a stronger repulsion from the interface than the repulsion induced by the intra-lattice potential energy barrier far from critical , thus enforcing the perturbed strongly pinned soliton component solitons to shed off a part of their power. A consequence is the formation of more stable interface asymmetric breathing modes with a smaller power.

The difference between the interface and intra-lattice potential energies may also explain an enhanced stability of symmetric complexes which are centered farther away from the interface. An example is plotted in Fig. 10, where purely real EVs are presented for symmetric soliton complexes whose components are centered at distances from the lattice interface in the corresponding lattices. Eventually, the symmetric complexes centered at are stable in their entire existence region. We stress that only the soliton complexes of the symmetric type were found, trying to set the centers of the component solitons farther from the interface.

The EV spectra for the asymmetric isosceles solitons, which are labeled by subscripts and in Figs. 3 and 6, indicate the exponential instability in their entire existence regions, see Fig. 8(b). This conclusion is confirmed by direct simulations, which show that perturbed asymmetric complexes of these types radiate a significant part of their power in the form of a breathing complex which moves across the corresponding chains, see Fig. 9(c). The instability of asymmetric solutions can be associated with the trend of the system to relax toward an energetically preferable state in the presence of the two above-mentioned forces – the repulsion from the lattice interface, and the force induced by the bulk lattice, which is measured by the respective Peierls-Nabarro potential barrier nashi (). The antisymmetric solutions with are shown to be unstable against oscillatory perturbations in the region of , see Fig. 8(c), while a narrow region with exponentially unstable antisymmetric complexes is found near the upper boundary for their existence domain. Dynamical calculations confirm the predictions of the linear stability analysis.

The isosceles soliton branch labeled by subscript in Figs. 3 and 6, which is stable in the whole existence region according to the linear-stability analysis, is stable in direct simulations as well. Therefore, it can be concluded that the symmetry-breaking bifurcation at is related to the stability exchange between the destabilized symmetric complex and the two emerging isosceles asymmetric complexes, one of which is stable and the other exponentially unstable. In addition, the antisymmetric complexes change their stability in the neighborhood of the bifurcation point (from the oscillatory instability at to the stability at ).

Returning to the global existence diagrams, it is worthy to note that two areas with coexisting symmetric, asymmetric and antisymmetric solutions can be identified: the domain featuring the coexistence of stable symmetric, unstable asymmetric and stable antisymmetric complexes, and the domain where the symmetric, antisymmetric and two asymmetric complexes are unstable and one isosceles solution branch is stable.

As said above, the dynamics of the interface soliton complexes is strongly related to the balance of the interface and bulk-lattice potential energies. When the inter-lattice coupling is too small, the interface complexes of both the symmetric and asymmetric types formed by the fundamental solitons are unstable. The exception is (as actually mentioned above) the isosceles asymmetric complex with (mode in Fig. 3 and 6), which is stable in its entire existence region. The strong repulsion from the interface, in comparison to the force from the bulk lattice, makes the formation of stable surface solitons difficult. By increasing the inter-lattice coupling, the energy difference between the interface and intra-lattice forces gets smaller, which makes it possible to create stable localized interface symmetric and antisymmetric complexes. With the further increase of the inter-lattice linkage against the intra-lattice coupling, which means proceeding to for fixed , the highly unstable asymmetric complex and stable antisymmetric one are the only soliton species generated by the system. These findings are correlated with results reported for the two-chain version of the present system in Ref. symnash (), where the stable symmetric, antisymmetric, and antisymmetric branches have been found in a certain part of the corresponding existence regions.

## 5 Conclusion

In this paper, we have analyzed the properties of fundamental localized interface complexes excited at the interconnection of three nonlinear semi-infinite chains. This system may be realized in nonlinear optics and BEC. In the framework of the system of three semi-infinite DNLS equations with the self-focusing on-site nonlinearity, coupled at the single site, we have found, by means of the VA (variational approximation) and numerical calculations, the threshold value of the inter-lattice coupling constant bounding the creation of symmetric surface soliton complexes. The existence of families of symmetric and asymmetric interface soliton complexes has been demonstrated. The variational predictions for the stationary modes were checked numerically. The stability analysis shows that the symmetric soliton complexes, which are created as a stable solution branch at the critical value of the lattice inter-coupling parameter, destabilize at the bifurcation point, where two isosceles asymmetric soliton complexes are created, one of them stable and the other exponentially unstable. In addition, the stability of the antisymmetric complexes changes twice – at the bifurcation point and near the upper boundary of their existence region. The third isosceles solution branch is exponentially unstable in its entire existence region. In other words, in the trilete system, the symmetric solution branch, three asymmetric branches and the antisymmetric one coexist in a part of the parameter space (past the symmetry-breaking bifurcation). These solution branches exist with their mirror-image counterparts. Direct simulations demonstrate that unstable symmetric complexes are transformed, radiating away a part of their power, into robust oscillating modes in the form of localized interface breathing complexes, while the unstable asymmetric complexes form breathers traveling away from the interface through the lattice. The origin of this behavior of the surface modes is related to the balance between the repulsion from the surface and the intra-lattice interactions.

Because the instability of asymmetric soliton complexes is stipulated by the repulsion from the interface (the linkage site, ), it seems plausible that they may be stabilized by making the on-site nonlinearity stronger at this site. Another interesting possibility is to consider solitons of the vortex type (cf. vortices studied in other linearly-coupled tri-core systems 20 (); Skryabin ()). The respective ansatz may be taken as , cf. ansatz (4). This vortex may be classified as one of the ”off-site” type, as its virtual pivot is located between the sites.

## Acknowledgments

Lj.H. and A.M. acknowledge support from the Ministry of Education and Science, Serbia (Project III45010).

## References

## References

- (1) I. E. Tamm, Z. Phys. 76 (1932) 849.
- (2) K. G. Makris, S. Suntsov, D. N. Christodoulides, G. I. Stegeman, and A. Hache, Opt. Lett. 30 (2005) 2466.
- (3) S. Suntsov, K. G. Markis, D. N. Christodoulides, G. I. Stegeman, A. Hache, R. Morandotti, H. Yang, G. Salamo, M. Sorel, Phys. Rev. Lett. 96 (2006) 063901.
- (4) F. W. Ye, Y. V. Kartashov, and L. Torner, Phys. Rev. A 74 (2006) 063616; F. K. Abdullaev, R. M. Galimzyanov, M. Brtka, L. Tomio, ibid. A 79 (2009) 056220; S. L. Cornish, N. G. Parker, A. M. Martin, T. E. Judd, R. G. Scott, T. M. Fromhold, C. S. Adams, Physica D 238 (2009) 1299; Y. V. Bludov, Z. Y. Yan, V. V. Konotop, Phys. Rev. A 81 (2010) 063610.
- (5) P. G. Kevrekidis, editor, The Discrete Nonlinear Schrödinger Equation: Mathematical Analysis, Numerical Computations, and Physical Perspectives, Springer: Berlin and Heidelberg, 2009.
- (6) M. Molina, R. Vicencio, Yu. S. Kivshar, Opt. Lett. 31 (2006) 1693.
- (7) Y. V. Kartashov, L. Torner, Opt. Lett. 31 (2006) 2172.
- (8) I. L. Garanovich, A. A. Sukhorukov, Yu. S. Kivshar, M. Molina, Opt. Express 14 (2006) 4780.
- (9) X. D. Cao, B. A. Malomed, Phys. Lett. A 206 (1995) 177.
- (10) W. Królikowski, Yu. S. Kivshar, J. Opt. Soc. Am. B 13 (1996) 876.
- (11) Q. E. Hoq, R. Carretero-González, P. G. Kevrekidis, B. A. Malomed, D. J. Frantzeskakis, Yu. V. Bludov, V. V. Konotop, Phys. Rev. E 78 (2008) 036605.
- (12) U. Peschel, R. Morandotti, J. S. Aitchison, H. S. Eisenberg, Y. Silberberg, Appl. Phys. Lett. 75 (1999) 1348.
- (13) R. Morandotti, H. S. Eisenberg, D. Mandelik, Y. Silberberg, D. Modotto, M. Sorel, C. R. Stanley, J. S. Aitchison, Opt. Lett. 28 (2003) 834.
- (14) M. Heinrich, R. Keil, F. Dreisow, A. Tünnermann, S. Nolte, A. Szameit, New Journal of Physics 12 (2010) 113020.
- (15) Y. J. He, W. H. Chen, H. Z. Wang, B. A. Malomed, Opt. Lett. 32 (2007) 1390.
- (16) E. Smirnov, M. Stepić, C. E. Ruter, D. Kip, V. Shandarov, Opt. Lett. 31 (2006) 2338.
- (17) D. N. Christodoulides, R.I. Joseph, Opt. Lett. 13 (1988) 794.
- (18) G. A. Siviloglou, K. G. Makris, R. Iwanow, R. Schiek, D. N. Christodoulides, G. I. Stegeman, Y. Min, W. Sohler, Opt. Exp. 14 (2006) 5508.
- (19) Y. V. Kartashov,V. A. Vysloukh, L. Torner, Phys. Rev. Lett. 96 (2006) 073901; C. R. Rosberg, D. N. Neshev, W. Królikowski, A. Mitchell, R. A. Vicencio, M. I. Molina, Yu. S. Kivshar, Phys. Rev. Lett. 97 (2006) 083901.
- (20) M. I. Molina, I. L. Garanovich, A. A. Sukhorukov, Y. S. Kivshar, Opt. Lett. 31 (2006) 2332.
- (21) A. Gubeskys, B. A. Malomed, I. M. Merhasin, Phys. Rev. A 73 (2006) 023607; S. K. Adhikari, B. A. Malomed, ibid. 77 (2008) 023607.
- (22) Y.-J. He, B. A. Malomed, in: P. G. Kevrekidis, editor, The Discrete Nonlinear Schrödinger Equation: Mathematical Analysis, Numerical Computations, and Physical Perspectives, Springer: Berlin and Heidelberg, 2009.
- (23) M. I. Molina, Yu. S. Kivshar, Opt. Lett. 33 (2008) 917.
- (24) V. A. Brazhnyi, V. V. Konotop, Mod. Phys. Lett. B 18 (2004) 627.
- (25) J. C. Eilbeck, P. S. Lomdahl, A. C. Scott, Physica D 16 (1985) 318.
- (26) E. M. Wright, G. I. Stegeman, S. Wabnitz, Phys. Rev. A 40 (1989) 4455.
- (27) A. W. Snyder, D. J. Mitchell, L. Poladian, D. R. Rowland, Y. Chen, J. Opt. Soc. Am. B 8 (1991) 2102.
- (28) C. Paré, M. Fłorjańczyk, Phys. Rev. A 41 (1990) 6287 ; A. I. Maimistov, Kvant. Elektron. 18 (1991) 758 [Sov. J. Quantum Electron. 21 (1991) 687]; B. A. Malomed, I. Skinner, P. L. Chu, G. D. Peng, Phys. Rev. E 53 (1996) 4084.
- (29) N. Akhmediev, A. Ankiewicz, Phys. Rev. Lett. 70 (1993) 2395; P. L. Chu, B. A. Malomed, G. D. Peng,J. Opt. Soc. A 10 (1993) 1379; J. M. Soto-Crespo, N. Akhmediev, Phys. Rev. E 48 (1993) 4710.
- (30) W. C. K. Mak, B. A. Malomed, P. L. Chu, J. Opt. Soc. Am. B 15 (1998) 1685.
- (31) A. Gubeskys, B. A. Malomed, Eur. Phys. J. D 28 (2004) 283.
- (32) A. Gubeskys, B. A. Malomed, Physical Review A 75 (2007) 063602; M. Matuszewski, B. A. Malomed, M. Trippenbach, ibid. 75 (2007) 063621; A. Gubeskys, B. A. Malomed, ibid. 76 (2007) 043623; M. Trippenbach, E. Infeld, J. Gocalek, M. Matuszewski, M. Oberthaler, B. A. Malomed, Phys. Rev. A 78 (2008) 013603; N. V. Hung, P. Ziń, M. Trippenbach, B. A. Malomed, Phys. Rev. E 82 (2010) 046602.
- (33) W. C. K. Mak, B.A. Malomed, P. L. Chu, Phys. Rev. E 55 (1997) 6134; ibid. E 57 (1998) 1092.
- (34) R. Driben, B. A Malomed, P. L. Chu, J. Phys. B: At. Mol. Opt. Phys. 39 (2006) 2455; L. Albuch, B. A. Malomed, Math. Comp. Simul. 74 (2007) 312; Z. Birnbaum, B. A. Malomed, Physica D 237 (2008) 3252; N. Dror, B. A. Malomed, ibid. 240 (2011) 526.
- (35) B. A. Malomed, Progr. Optics 43 p. 71 (E. Wolf, editor: North Holland, Amsterdam, 2002).
- (36) G. Herring, P. G. Kevrekidis, B. A. Malomed, R. Carretero-González, D. J. Frantzeskakis, Phys. Rev. E 76 (2007) 066606.
- (37) C. Chong, D. E. Pelinovsky, Discrete and Continuous Dynamical Systems S 4 (5) (2011) 1019, doi:10.3934/dcdss.2011.4.1019.
- (38) Lj. Hadžievski, G. Gligorić, A. Maluckov, B. A. Malomed, Phys. Rev. A 82 (2010) 033806.
- (39) A. Szameit, Y. V. Kartashov, D. Dreisow, T. Pertsch, S. Nolte, A. Tünnermann, L. Torner, Phys. Rev. Lett. 98 (2007) 173903; A. Szameit, Y. V. Kartashov, D. Dreisow, M. Heinrich, V. A. Vysloukh, T. Pertsch, S. Nolte, A. Tünnermann, F. Lederer, L. Torner, Opt. Lett. 33 (2008) 663; A. Szameit, Y. V. Kartashov, V. A. Vysloukh, M. Heinrich, F. Dreisow, T. Pertsch, S. Nolte, A. Tünnermann, F. Lederer, L. Torner, ibid. 33 (2008) 1542; A. Szameit, I. L. Garanovich, M. Heinrich, A. A. Sukhorukov, F. Dreisow, T. Pertsch, S. Nolte, A. Tünnermann, Y. S. Kivshar, Phys. Rev. Lett. 101 (2008) 03902; M. Heinrich, Y. V. Kartashov, L. P. R. Ramirez, A. Szameit, F. Dreisow, R. Keil, S. Nolte, A. Tünnermann, V. A. Vysloukh, L. Torner, Opt. Lett. 34 (2009) 3701.
- (40) A. Sigler, B. A. Malomed, D. V. Skryabin, Phys. Rev. E 74 (2006) 066604.
- (41) B. A. Malomed, M. I. Weinstein, Phys. Lett. A 220 (1996) 91; D. J. Kaup, Mathematics and Computers in Simulations 69 (2005) 322; R. Carretero-González, J. D. Talley, C. Chong, B. A. Malomed, Physica D 216 (2006) 77.
- (42) Lj. Hadžievski, G. Gligorić, A. Maluckov, B. Malomed, Phys. Rev. A 82 (2010) 033806.
- (43) D. J. Kaup, B. A. Malomed, Physica D 184 (2003) 153.
- (44) Y. Sivan, B. Ilan, and G. Fibich, Phys. Rev. E 78 (2008) 046602.
- (45) G. Gligorić, A. Maluckov, Lj. Hadžievski, B. A. Malomed, Phys. Rev. A 79 (2009) 053609; G. Gligorić, A. Maluckov, Lj. Hadžievski, B. A. Malomed, Phys. Rev. A 78 (2008) 063615.