# Peierls-Nabarro energy surfaces and directional mobility of discrete solitons in two-dimensional saturable nonlinear Schrödinger lattices

###### Abstract

We address the problem of directional mobility of discrete solitons in two-dimensional rectangular lattices, in the framework of a discrete nonlinear Schrödinger model with saturable on-site nonlinearity. A numerical constrained Newton-Raphson method is used to calculate two-dimensional Peierls-Nabarro energy surfaces, which describe a pseudopotential landscape for the slow mobility of coherent localized excitations, corresponding to continuous phase-space trajectories passing close to stationary modes. Investigating the two-parameter space of the model through independent variations of the nonlinearity constant and the power, we show how parameter regimes and directions of good mobility are connected to existence of smooth surfaces connecting the stationary states. In particular, directions where solutions can move with minimum radiation can be predicted from flatter parts of the surfaces. For such mobile solutions, slight perturbations in the transverse direction yield additional transverse oscillations with frequencies determined by the curvature of the energy surfaces, and with amplitudes that for certain velocities may grow rapidly. We also describe how the mobility properties and surface topologies are affected by inclusion of weak lattice anisotropy.

###### pacs:

42.65.Wi, 63.20.Pw, 63.20.Ry, 05.45.Yv## I Introduction

Intrinsically localized modes, or discrete solitons (breathers), appear as generic excitations in a large variety of physical systems rep2 (), where spatial periodicity (or discreteness) provides gaps in the linear dispersion relation and nonlinearity allows for detuning the oscillation frequencies into these gaps, resulting in spatial localization. Particularly important applications of discrete nonlinear systems are nonlinear optical waveguide arrays rep1 (), with possibilities to change and control all essential parameters, such as geometry, dimensionality, nonlinearity, beam angle, etc. For weakly coupled waveguides, the relevant mathematical description, derived via coupled mode theory cj88 (); rep1 (), is a Hamiltonian lattice equation of the discrete nonlinear Schrödinger (DNLS) type rep2 (); rep1 (); cj88 (); ej03 (), where the particular type of nonlinearity depends on the nonlinear response of the medium. The most studied case is a DNLS equation with cubic on-site nonlinearityej03 (), corresponding to a Kerr medium rep1 (), which also appears generically as a modulational equation for the small-amplitude dynamics in chains of coupled anharmonic oscillators rep2 (). However, photorefractive media also enables the generation of discrete spatial solitons rep1 (), and the corresponding lattice model is a DNLS equation with saturable nonlinearity (s-DNLS) which was studied in a number of theoretical works prlkip (); jpakhare (); plaeil (); prlmel (); ob07 (); pre73 (); rnkf09 ().

A particularly interesting property of the one-dimensional (1D) s-DNLS model prlkip (); jpakhare (); plaeil (); prlmel (); ob07 () is the existence of certain “sliding velocities”, where localized discrete solitons may travel in the lattice without radiation. This behavior was connected to the existence of “transparent points” associated with the vanishing of a so-called Peierls-Nabarro (PN) potential barrier pre48 (), usually defined as the difference in energy (Hamiltonian) at constant power (norm) between the two fundamental localized stationary solutions centered at one site (odd mode) and symmetrically in between two sites (even mode). Close to the points of vanishing of this energy difference are regions of stability exchange between the even and odd solutions, associated in general with the appearance of a family of intermediate, asymmetric stationary solutions, connecting both types of symmetric solutions at the bifurcation points pdaub (). The existence of regimes of enhanced mobility close to such bifurcation points is also well-known from studies of other one-dimensional lattice models pdaub (); pdcret (); oje03 ().

On the other hand, it is still an open issue whether moving discrete solitons may exist as localized, radiationless modes also in two-dimensional (2D) lattices. As was shown numerically in pre73 (), a scenario with exchange of stability through bifurcations with asymmetric stationary solutions appears also for the 2D saturable model in a square (isotropic) lattice, involving in this case three different types of fundamental solutions krb00 (): one-site (odd-odd, OO), two-site (odd-even, OE), and four-site (even-even, EE) modes. It was also shown numerically in pre73 (), that solutions with good (but generally not radiationless) mobility in the axial directions may exist in these regimes, and that the necessary energy needed for rendering a mobile stable stationary solution agreed well with the concept of PN-barrier, if its definition was extended to take into account also the energy for the relevant intermediate stationary solution (an analogous situation is well-known for the PN potential of kinks pr82 ()). Very similar observations were also made later for a 2D DNLS model with cubic-quintic nonlinearity ccmk09 (). Moreover, in the low-power regime of the 2D s-DNLS equations, good mobility was also observed in diagonal directions pre73 ().

However, the relation between existence of regions of stability exchange, small PN-barrier and mobile localized solutions in 2D is not that trivial, as shown, e.g., for a model with cubic inter-site nonlinearities in oj09 (): even in regimes with small PN barrier and existence of intermediate solutions, the mobility may be very poor, if there is no continuous path in phase space passing close to the relevant stationary modes. On the other hand, as was observed for the low-power (i.e., close to continuum limit) regime of a 2D lattice with quadratic nonlinearity in susanto07 (), the effective Peierls-Nabarro potential may in some situations be weak enough to allow mobility in arbitrary directions, without any direct connection to bifurcations and symmetry-broken stationary solutions.

So there is clearly need for a better understanding of the conditions for mobility in 2D lattices. It is the purpose of the present paper to generalize the concept of PN-barrier as discussed above, and introduce a full 2D PN potential surface describing the pseudo-potential landscape in-between all stationary modes. We will use a numerical constrained Newton-Raphson (NR) method, previously applied to 1D lattices in rodrigo1d (), to explicitly construct these surfaces for the 2D saturable model from pre73 (), and show how parameter regimes and directions of good mobility may be immediately identified from smooth, flat parts of these surfaces. We will also illustrate how the interplay between translational motion in one lattice direction and oscillatory motion in the orthogonal direction can be intuitively understood from the topology of the corresponding PN surfaces.

The structure of the paper is as follows. In Sec. II, we describe the 2D s-DNLS model, its stationary localized solutions (Sec. II.1), and the explicit implementation of the numerical constraint method used to calculate the PN potential surfaces (Sec. II.2). We point out some essential differences between the 2D implementation and earlier implementations for the 1D case. In Sec. III.1 we identify the regimes of parameter space where smooth surfaces can be found for the isotropic lattice, and make a classification of appearing surfaces with different topologies. In Sec. III.2 we illustrate with direct dynamical simulations various types of behavior for discrete solitons moving slowly in different directions, corresponding to the surfaces discussed in Sec. III.1. We also here show examples on how simultaneous excitation of transverse oscillation modes may affect discrete solitons moving slowly in axial directions. In Sec. IV, we analyze effects of adding a weak anisotropy to the lattices in the parameter regimes used in previous sections, in order mainly to investigate whether this may promote mobility in directions different from the axial or diagonal ones (which presumably would be associated with the appearance of additional ”valleys” in the PN potential surface). Finally, in Sec. V we make some concluding remarks, summarizing our results and pointing out some relevant directions for future research.

## Ii Model

We consider the following general form of the 2D s-DNLS equation in a rectangular geometry pre73 ():

(1) |

where corresponds to the normalized propagation coordinate, to an effective nonlinear parameter, and represents the (complex) electric field amplitude at site in an array of sites. defines the linear interaction between nearest-neighbor waveguides, where is an anisotropy parameter allowing us to study different 2D scenarios from a completely isotropic lattice () to a completely decoupled set of 1D arrays (). Model (1) possesses two conserved quantities, the power (norm) defined as

(2) |

and the Hamiltonian (energy) defined as

(3) |

### ii.1 Stationary solutions

Stationary solutions to Eq. (1) are of the form , where is a -independent (generally complex) amplitude and corresponds to the propagation constant or frequency pre73 (). Extended stationary solutions with constant amplitude (plane waves) exist in frequency bands, whose edges depend on the amplitude . The linear band region is obtained for low-amplitude plane waves (), and it is easy to show that they will exist in the region . Moreover, for high-amplitude plane waves (), the nonlinear term saturates and the corresponding interval is .

The fundamental localized stationary solutions to Eq. (1), defined as solutions with real amplitudes having a single maximum distributed on the sites in one unit-cell of the lattice, were described for the isotropic case in pre73 (). In general, we may identify one-site (OO), two-site horizontal (EO), two-site vertical (OE), and four-site (EE) solutions. Figs. 1(a)-(d) show typical profiles for these types of excitations. As discussed in pre73 (), these in-phase localized solutions bifurcate from the fundamental linear mode (upper band edge) in the small-amplitude limit, and merge into the corresponding high-amplitude plane-wave mode in the infinite-power limit. As a consequence, they exist in the region , so that the size of the existence region will be just . [A two-site diagonal solution Kalosakas () [see Fig. 1(h)] could also be excited but only in frequency regimes where solutions are very localized with , requiring large values of (for , )]. In addition, in certain parameter regimes there are also asymmetric, intermediate stationary solutions, associated with exchange of stability between the symmetric ones pre73 (). We may identify three types of such solutions, termed henceforth intermediate 1 (IS1), intermediate 2 (IS2), and intermediate 3 (IS3), respectively. IS1 and IS3 both connect a stable two-site mode, OE/EO, with another simultaneously stable solution: the one-site mode OO (IS1), or the four-site mode EE (IS3) [see Figs. 1(e) and (g)]. Thus, in these cases the unstable intermediate solutions act as carriers of instability between the corresponding fundamental modes. The unstable IS2 solution exists when the two-site solutions are stable. It connects the unstable one-site solution OO with the likewise unstable four-site solution EE, and stabilizes the latter mode when reaching it [see Fig. 1(f)].

Figs. 2 shows fundamental properties of different stationary solutions for an isotropic lattice with (similar curves were plotted and discussed in pre73 (), but we here include some additional examples to facilitate the comparison with the full PN surfaces in Sec. III.1). Power versus frequency curve is shown in Fig. 2 (a), where different solutions repeatedly cross each other in the whole range of parameters. In Fig. 2 (a)-inset we plot as a function of , for . This figure shows clearer how solutions cross each other including many intermediate solutions appearing in the region . Moreover, we also observe these crossings when plotting versus power [Fig. 2 (b)]. (We plot instead of because Hamiltonian differences between stationary solutions of the same power are generally quite small, which is indeed a favorable scenario for moving solutions). Inset in Fig. 2 (b) shows a zoom in the region where the Hamiltonian of different fundamental solutions matches at some points. This was initially interpreted as a vanishing PN barrier prlkip (), which is actually not the case due to the intermediate solutions appearing in-between, raising the effective energy barrier pre73 () [see Fig. 2 (b)-inset when thick blue lines coincide]. However, good mobility would be expected close to these regions but, still, it will be strongly determined by the specific kick or perturbation given to the solution in order to put it in movement.

In order to study the stability of stationary localized solutions we implement a standard method jpakhare (); pre73 (), obtaining all the stable/unstable linear perturbation-modes and their corresponding instability measure , where are the oscillation frequencies of the linear eigenmodes. Thus, if the solution is (at least marginally) stable and if it is unstable. Stability versus power for the same solutions as in Figs. 2 (a), (b) is shown in Fig. 2 (c). (Note that in the corresponding figures 2-3 of pre73 (), the quantity was used as instability measure.) From this figure, we can see how unstable intermediate solutions IS1, IS2, IS3 appear when two or three solutions (regarding OE and EO as different solutions) are simultaneously stable. Note that we never observed regions for simultaneously stable one and four-site solutions for isotropic coupling. However, as will be discussed further in Sec. IV, in the anisotropic case such regions exist.

### ii.2 Constraint method

The constraint method allows us to construct energy surfaces connecting stationary solutions for a given value of power. In that sense, it helps us to effectively predict and interpret the dynamics across the lattice. Critical points will represent stationary solutions, and a coherent movement across the lattice should transform one solution to the other by keeping the power constant pre48 (). The method was originally introduced by Aubry and Cretegny pdaub () and then implemented by Savin et al. savin () to study the mobility of kinks in nonlinear Klein-Gordon lattices (see also Ref. Sepulchre () for a related approach to analyze travelling breathers in 1D oscillator chains in terms of an effective Hamiltonian.) Lately, it was numerically implemented to analyze surface states in one-dimensional semi-infinite systems rodrigo1d (). By adiabatically changing the amplitude in one particular site, specifically chosen as the one after the main peak (the peak is at and the constrained amplitude at ), the one- and the two-site solutions could be connected. Ending the sweep when and the center of mass of the constrained solution, , has varied from to , a one-dimensional energy surface, vs , can be sketched. Technically speaking, the method used in Ref. rodrigo1d () consists on eliminating one equation from the Newton-Raphson problem, the one of the constrained amplitude which is not anymore an unknown variable. However, as the power is kept constant, an equation for is added and the frequency becomes a variable completing the variable-equations set.

In the present work, for the 2D lattice we implement a more sophisticated method where we explicitly vary the center of mass instead of the amplitude. Since we are interested in the energy landscape around the fundamental stationary solutions, we will assume also the amplitude of the constrained solutions to be real and positive on the constraint sites. Then, from the definition of the center of mass coordinates in the horizontal and vertical directions,

(4) |

we can solve for any pair of amplitudes and at places and :

(5) | |||||

where with . Thus, as before the actual constraints will be in the amplitudes and , but now we can tune the center of mass as wished, from a given stationary solution towards any other. In general, the constrained solutions obtained in this way will of course not be stationary solutions of the full system. To identify any stationary solution - including the intermediate ones - we check whether the value of obtained from the constrained NR scheme coincides with the frequency of a hypothetical true stationary solution to the full Eq. (1) with the computed amplitude profile. Furthermore, constrained solutions allow us to calculate their Hamiltonian, and therefore to construct an effective energy landscape.

The choice for positions and is evidently not unique, but in order to most efficiently trace out a smooth energy landscape connecting the fundamental solutions (if it exists), the constraint sites should preferably be chosen within the unit cell where the amplitudes are large. It turned out that, starting from a stationary one-site solution centered at , in most cases the best option is to choose the first constraint at site or , and the other one at site , as sketched in Fig. 3. The results shown in the following sections are obtained using these constraint sites. We also tried using constraint sites at and (i.e., along a diagonal); however, with this choice we typically were not able to find the four-site EE solution when starting from the one-site OO. Instead, the NR method with constraints on these sites generally converges, at , to a two-site diagonal solution (which as remarked in Sec. II.1 is an exact solution to the unconstrained equation only for large ), which is not expected to be relevant for the mobility process.

## Iii Isotropic lattice

### iii.1 Description of PN surfaces

By exploring the parameter-space , we have identified regions where we were able to compute complete two-dimensional energy surfaces (using constraint as shown in Fig. 3) and regions where we cannot. In general, complete surfaces could not be obtained for large values of the nonlinearity constant, , when solutions get strongly localized. In some cases we could trace some specific mobility directions (e.g., along lattice axes or diagonal) connecting stationary solutions with almost equal Hamiltonian, but not the full two-dimensional scenario involving all fundamental solutions. Note also that, as remarked in pre73 (), for large and small power the s-DNLS model essentially behaves as a cubic DNLS model with effective nonlinearity . Consequently, we could not find complete energy surfaces with all fundamental solutions for the cubic DNLS either (although surfaces involving the two-site diagonal solution in place of the EE solution could be obtained as mentioned in Sec. II.2; for the cubic DNLS this two-site solution does exist as a true stationary solution above some threshold level of power Kalosakas ()).

On the other hand, for very small values of , the computation of energy surfaces becomes difficult for technical reasons: due to the widening of the solutions, considerably larger lattices are needed to remove the influence from boundary effects. Therefore, in the following, we will present the main phenomenology for the isotropic case () found for intermediate values of , , where complete surfaces involving all four fundamental solutions are obtained for all values of the power. We will show results for the particular value , but the scenario is found to be qualitatively the same for all in this interval.

From Fig. 2 (c), we may identify essentially five different regimes where energy surfaces of qualitatively different nature should be expected, depending on the level of power. The first one is for low power, where (similarly to the cubic DNLS) the one-site solution is always stable and the other fundamental stationary solutions are all unstable. The corresponding energy surface is illustrated in Fig. 4 (a), where the one-site solution yields the energy minimum, the two-site solutions saddle points and the four-site solution the maximum. Note that, in this low-power regime, the surfaces for this value of are still rather flat, and therefore some mobility may result if the one-site solution is kicked to overcome the barriers, in the axial as well as in the diagonal directions, as illustrated by Figs. 4 (a) and (c) in Ref. pre73 (). (Another example of mobility in this regime is discussed in Sec. III.2 below.)

The saturable nature of the system becomes evident at higher powers. In the second regime, appearing for the first time when , the one- and the two-site solutions are stable simultaneously and, as shown in Fig. 4 (b), these three points all correspond to local minima of the surface. Intermediate solutions IS1 connecting the one- with the two-site solutions in the horizontal and vertical directions appear as saddle points [white dots in Fig. 4 (b)]. Note that the energy landscape for the parameter values in Fig. 4 (b) is almost flat between the one-site and two-site solutions in the axial directions, leading to the very good axial mobility shown in Fig. 4 (d) of Ref. pre73 (), while the maximum corresponding to the unstable four-site solution creates a too large effective barrier to overcome in the diagonal direction.

The third power region is the one in which only the two-site solutions are stable. It appears for the first time when , and is illustrated in Fig. 4 (c). Here, the stable two-site solutions correspond to two local minima of the surface, and the unstable one- and four-site solutions both to local maxima. The two unstable solutions are connected by the intermediate solution IS2 (white dot at the surface), corresponding to a saddle point which for symmetry reasons (for the isotropic lattice) will lie along the diagonal connecting the unstable solutions. As will be illustrated in Sec. III.2, the easiest mobility in this case is expected to occur in a diagonal direction, connecting the two stable stationary solutions.

The fourth regime corresponds to simultaneously stable two- and four-site solutions. Such a regime appears for the first time when . If we further increase the power, we will find a similar region with these three simultaneously stable solutions for . In this regime, as illustrated in Fig. 4 (e), the two- and four-site solutions all correspond to minima, the unstable intermediate solutions IS3 to saddle points, and the one-site solution to a maximum. Thus, comparing Figs. 4 (b) and (e), we see that the structure of the potential has been completely inverted. Now, a big “hill” is located at the one-site solution, and as a consequence no mobility is expected involving this solution. As will be shown below in Sec. III.2, the simplest mobility scenario will now be the one in which the two-site solutions travel through the four-site one across the lattice.

The fifth regime, with the four-site solution being the only stable one, occurs for the first time when . As illustrated in Fig. 4 (d), the four-site solution is now a minimum of the PN potential, while the other three unstable solutions correspond to saddle points (EO and OE) and maximum (OO), respectively. (There are no intermediate solutions in this regime.) Thus, by increasing the power we have now reverted the surface compared to the low-power regime in Fig. 4 (a).

Further increasing of power shows for a new region of stable two-site solutions, corresponding to the third regime. For the one- and two-site solutions are simultaneously stable, so the scenario is equivalent to the second regime, followed again for by the first regime, respectively. This second complete inversion thus corresponds to a regain of the low-power characteristics of the surface. Furthermore, repetition of these scenarios can be found for , and we confirmed, e.g., the existence of a complete, smooth surface, analogous to the one shown in Fig. 4 (d) for the fifth regime, for .

### iii.2 Mobility dynamics for isotropic lattice

To explicitly show the connection between the different types of energy surfaces described in Sec. III.1 and the mobility of localized solutions, we numerically integrate model (1) by taking as initial condition a stationary solution perturbed with a small kick: , where and correspond to the kick strength in the horizontal and vertical directions, respectively. If the surfaces would be completely flat, solutions should move even with an infinitesimally small kick. However, as we have seen, in general surfaces are not flat due to the discreteness and the self-induced PN potential, and although the PN barrier can be very small in certain directions, it is generally non-zero. Therefore, in order to put a localized solution in movement, some amount of kinetic energy (represented by this kick) should be given to effectively overcome the energy barriers.

We first discuss an example from the first, low-power, regime corresponding to the energy surface in Fig. 4 (a) (as mentioned above, different mobility examples from this regime were shown in Figs. 4 (a) and (c) in Ref. pre73 ()). Here, we took the (unstable) EO solution and kicked it with a very small value in both directions, yielding the dynamics numerically observed in Fig. 5 (a). In the -direction, the kinetic energy is not sufficient to overcome the barrier created by the four-site solution, while in the -direction it can move towards the minimum corresponding to the one-site solution. The resulting dynamics for the center-of-mass positions shows how the movement in the horizontal direction gets combined with oscillations in the vertical one. The dependence vs also clearly shows how the solution feels the potential in terms of its velocity: maximum velocities occur around integer lattices sites (where the potential is a minimum) while the minimum velocities occur close to middle points [see full line in Fig. 5 (a); compare also with the analogous scenarios for the kicked OO mode in Figs. 4 (a) and (c) of Ref. pre73 ()].

In the second power regime, with energy surfaces as in Fig. 4 (b), good mobility can be expected only in axial directions, with effective energy barrier determined by the intermediate solution between the one- and the two-site solutions as discussed in Ref. pre73 (). The dynamics in this power regime, with a very small kick applied only in the axial direction to one of the stable stationary solutions, will be as already illustrated in Fig. 4 (d) of Ref. pre73 (): the solution moves very slowly and adiabatically traces the shape of the potential with a minimal velocity at the places of the intermediate solutions.

In the third power regime, with surfaces as in Fig. 4 (c), a very interesting kind of mobility is observed: a diagonal mobility between the (stable) horizontal and vertical two-site solutions as illustrated in Fig. 5 (b). The initial EO solution, kicked equally in the - and -directions, gets sufficient kinetic energy to pass over the small barrier created by the intermediate IS2 solution. It then continues through the OE solution, passes another IS2 barrier, and then to the other EO solution shifted by one site in both directions. Although the potential connecting these two solutions is not completely flat, there is a very good transport of energy in this direction, allowing mobility for more than one lattice diagonal in the considered lattice.

If we take a look at the surface in the fifth power regime [Fig. 4 (d)], we realize that an initial (unstable) one-site solution may move in any direction by slightly kicking it since it corresponds to a maximum. Fig. 5 (c) shows an example for kicking the one-site solution symmetrically in both directions in order to make it move passing through the four-site, i.e., a diagonal movement. The velocity has maximum in the minima of the potential (corresponding to the EE solution) and minimum in the potential maxima (OO solution) (no intermediate solutions appear in this regime).

Finally, in the fourth power regime, we see from the surface [Fig. 4 (e)] that the two- and four-site solutions are stable simultaneously, presenting an intermediate solution in-between that will define the effective energy barrier. This barrier is very small and, therefore, a very small kick is also needed. Figure 5 (d) shows the evolution starting from a two-site vertical solution, passing through the intermediate one and arriving to the four-site solution.

As illustrated by the example in Fig. 5(a), the energy surfaces also provide nice intuitive interpretations to the dynamics of moving discrete solitons with additional perturbations transverse to the direction of motion. In this example, the dynamics in the two orthogonal directions appear essentially independent of each other (propagation in while oscillating in ). However, as we will illustrate with another example below, there are situations where the particular surface topologies may lead to a more intricate interplay between the oscillatory and translational dynamics.

We present in Fig. 6 a case with . With the notation from Sec. III.1, this value of the power belongs to the second power regime, and the structure of the surface is phenomenologically identical to the one sketched in Fig. 4(b), but for a higher level of power. A picture of this energy landscape, periodically extended along the -direction, is shown in Fig. 6 (a). In this figure, the darkest regions correspond to the positions for the minima (the OO and EO/OE modes), while the brightest regions correspond to maxima (EE modes). Note also the positions of the saddle points (IS1, marked with arrows) in-between the OO and EO/OE modes.

The kick we choose for the vertical direction is very small, , in order to disturb the horizontal movement only weakly. Since the potential shape in the -direction - in a first approximation - will correspond to a harmonic potential, the profile will essentially perform harmonic oscillations in the -direction, with a frequency (determined by the curvature of the surface) that will vary only slightly as the mode translates in the -direction. To look for possible resonances between the oscillatory and translational dynamics, we implement the following scheme: For fixed and we numerically integrate model (1), starting with the OO stationary solution centered at position , and measure the z-values for which is equal to and , respectively. For these z-values we compute , and plot in Fig. 6(b) as a function of the horizontal kick strength .

The particular values of where symbols in Fig. 6(b) coincide for correspond to kick strengths where the natural frequency of the surface potential is commensurate with the translational velocity. In Fig. 6(c) we show some examples for these particular points. For the solution makes two cycles before arriving to the next integer position in . For the solution makes one and a half cycles, for one cycle, for one half of cycle and for one quarter of cycle, before arriving to the next integer position in . As can be seen, around these points the oscillations apparently remain bounded and the dynamics stable, at least for long times. (Note that we do not see any visible phase-locking effects between the translational and oscillatory motion here: the symbols in Fig. 6(b) seem to coincide only at isolated points and not in intervals. This may be an indication that the dissipative effects on the effective dynamics are very weak.)

However, in other regimes we observe that the amplitude of the oscillations in the -direction rapidly grows. Two examples are illustrated in Fig. 6(d), corresponding to values of where the -displacement in Fig. 6(b) has maxima. (Here, we have used a larger system of sites, in order to clearly identify the nature of this particular dynamics.) The initial increase of the -amplitude suggests a kind of unstable dynamics. However, as is clear from the longer simulation in Fig. 6(d), the oscillation amplitude in the Y-direction remains bounded and thus the trajectory does not escape in the vertical direction but remains trapped by the barrier. The typical dynamics, thus, exhibits amplitude-modulated oscillations (beatings), with the largest amplitudes corresponding to maxima close to integer where the surfaces are flatter. Thus, this kind of mobility is a signature of the 2D topology of the potential and it also validates the surfaces computed with our method. (We performed several computations in different -regimes finding similar results as the ones presented in Fig. 6).

## Iv Effects of weak anisotropy

Considering weak anisotropy by putting in Eq. (1) (implying a stronger coupling in the horizontal than in the vertical direction), the main qualitative modifications as concerns the stationary solutions are shifts of the stability regions for the different solution types (evidently, the two-site horizontal, EO, and vertical, OE, solutions are now non-equivalent). Focusing our discussion on the parameter value analyzed in detail for the isotropic case in the previous sections, the results for the stability analysis are presented in Fig. 7. We find that, for fixed and increasing power, the first stability inversion appears between the one-site OO and the horizontal two-site EO solutions, and it occurs for lower values of the power than in the isotropic case. The first interval of stability of the four-site EE solution gets narrower when decreases, mainly because its lower limit increases (with a corresponding increase of the upper stability limit for the EO solution). The first stability region of the vertical two-site OE solution disappears for . On the other hand, its second stability regime for higher powers is enhanced, whereas the corresponding stability regime for its horizontal counterpart gets narrower and disappears for .

With stronger anisotropy, some qualitatively new regimes not present in the isotropic case will appear. For , a new region of stability of the OO mode appears around , leading to a small window of multistability where the OO, EO and EE modes are stable simultaneously. An example of an energy surface from this regime is shown in Fig. 8 (a), where also two intermediate solutions [white dots] can be found on the edges. There is also a regime with bistability of the OO and EE modes only, as illustrated by the pseudo-potential surface in Fig. 8 (b). Note that, in this case, the intermediate solution does not leave the edge, in contrast to the IS2 mode existing on the diagonal between the one- and four- site solutions for the isotropic case of simultaneously stable EO/OE solution [Fig. 4 (c)]. Since anisotropy destroys the symmetry between the OE and EO solutions, bistability of both 2-site solutions disappears already for a small amount of anisotropy, and for smaller the IS position can only be found on an edge of the potential surface.

Thus, since no additional “valleys” were found to develop in the PN surfaces (all minima appear on the edges), the best mobility for anisotropic lattices should generally be expected to occur along lattice directions (at least in cases where diagonal coupling terms can be neglected, as assumed throughout this work).

## V Conclusions

In conclusion, we have studied in a deeper way the problem of mobility of localized modes in two-dimensional saturable discrete systems. We numerically implemented a constrained Newton-Raphson method to construct full Peierls-Nabarro energy surfaces, which appeared as very useful tools for predicting the dynamical properties of localized excitations. Although these surfaces were never found to be completely flat (and therefore the corresponding Peierls-Nabarro barriers is strictly never zero), parameter regimes and directions of good mobility were seen to correspond to smooth, flat parts of the surfaces.

For the isotropic saturable model, five different surface topologies could be identified in different power regimes, depending on the stability properties of the different fundamental stationary solitons. By numerically studying the dynamics of perturbed stationary solutions, we showed how these different topologies yielded qualitatively different kinds of optimal mobility, generally in axial or diagonal directions and with velocities varying according to the shape of the potential while the profile is propagating across the lattice. The energy surfaces were also found to be very useful for interpreting the dynamics resulting from the interplay between translational and oscillatory motion in orthogonal directions.

We also studied the effect of weak anisotropy, where the breaking of the symmetry of the system also yields energy surfaces with lower symmetry. In addition, we showed how this symmetry breaking could lead to a change of stability properties for the fundamental stationary solutions, and as a result two new types of surface topologies were identified. From the shape of these surfaces, we predicted that the best mobility in the anisotropic case should generally occur along lattice directions.

Thus, the method developed here appears to be a very good tool for the understanding of the mobility dynamics of localized excitations in higher-dimensional lattices, clearly showing the “real” energy barriers that localized solutions experience. We hope that this approach will be useful also for other types of nonlinar lattice models.

Finally, from an application point of view, we have given further examples on how the saturable nature of nonlinearity promotes the generation of many different stability and mobility scenarios. The example presented in section III.2 shows the possibility to steer the dynamics of the solution by using the natural frequencies of its own self-induced potential. As an application of this, we propose a discriminative optical switch which, for example, will only be activated when and positions are - simultaneously - an integer number. This is one more example of nonlinear lattices being the key to success for controlling the propagation of light in future all-optical technologies.

###### Acknowledgements.

We are grateful to Serge Aubry and Sergej Flach who both, independently, suggested this approach to analyze 2D mobility, and also for discussions and comments on the manuscript. We thank Mario I. Molina for useful discussions. The authors acknowledge financial support from different sources: FONDECYT 1070897 and 7080001, FB0824/2008, and CONICYT fellowship. Authors acknowledge the hospitality of the Max Planck Institute for the Physics of Complex Systems (MPIPKS), Dresden, where this work was in part developed. M.J. also acknowledges support from the Swedish Research Council.## References

- (1) S. Flach and A.V. Gorbach, Phys. Rep. 467, 1 (2008).
- (2) F. Lederer et al., Phys. Rep. 463, 1 (2008).
- (3) D.N. Christodoulides and R.I. Joseph, Opt. Lett. 13, 794 (1988).
- (4) J.C. Eilbeck and M. Johansson, in Localization and Energy Transfer in Nonlinear Systems, Proceedings of the Third Conference, San Lorenzo de El Escorial Madrid, edited by L. Vázquez, R.S. MacKay, and M.P. Zorzano (World Scientific, Singapore, 2003), p.44.
- (5) M. Stepić, D. Kip, Lj. Hadžievski, and A. Maluckov, Phys. Rev. E 69, 066618 (2004); Lj. Hadžievski, A. Maluckov, M. Stepić, and D. Kip, Phys. Rev. Lett. 93, 033901 (2004); A. Maluckov, M. Stepić, D. Kip, and Lj. Hadžievski, Eur. Phys. J. B 45, 539 (2005); A. Maluckov, Lj. Hadžievski, and M. Stepić, Physica D 216, 95 (2006); Eur. Phys. J. B 53, 333 (2006).
- (6) A. Khare, K.Ø. Rasmussen, M.R. Samuelsen, and A. Saxena, J. Phys. A38, 807(2005); ibid. 42, 085002 (2009); ibid. 43, 375209 (2010).
- (7) J. Cuevas and J.C. Eilbeck, Phys. Lett. A 358, 15 (2006).
- (8) T.R.O. Melvin, A.R. Champneys, P.G. Kevrekidis, and J. Cuevas, Phys. Rev. Lett. 97, 124101 (2006); Physica D 237, 551 (2008).
- (9) O.F. Oxtoby and I.V. Barashenkov, Phys. Rev. E 76, 036603 (2007).
- (10) R.A. Vicencio and M. Johansson, Phys. Rev. E 73, 046602 (2006).
- (11) V.M. Rothos, H.E. Nistazakis, P.G. Kevrekidis, and D.J. Frantzeskakis, J. Phys. A 42, 025207 (2009).
- (12) Yu.S. Kivshar and D.K. Campbell, Phys. Rev. E 48, 3077 (1993).
- (13) T. Cretegny, Ph.D. thesis (in French), École Normale Supérieure de Lyon, 1998; S. Aubry, Physica D 216, 1 (2006).
- (14) S. Aubry and T. Cretegny, Physica D 119, 1 (1998).
- (15) M. Öster, M. Johansson, and A. Eriksson, Phys. Rev. E 67, 056606 (2003).
- (16) P.G. Kevrekidis, K.Ø. Rasmussen, and A.R. Bishop, Phys. Rev. E 61, R2006 (2000).
- (17) M. Peyrard and M. Remoissenet, Phys. Rev. B 26, 2886 (1982).
- (18) C. Chong, R. Carretero-González, B.A. Malomed, and P.G. Kevrekidis, Physica D 238, 126 (2009).
- (19) M. Öster and M. Johansson, Physica D 238, 88 (2009).
- (20) H. Susanto, P.G. Kevrekidis, R. Carretero-González, B.A. Malomed, and D.J. Frantzeskakis, Phys. Rev. Lett. 99, 214103 (2007).
- (21) M.I. Molina, R.A. Vicencio, and Yu.S. Kivshar, Opt. Lett. 31, 1693 (2006); C.R. Rosberg, D.N. Neshev, W. Krolikowski, A. Mitchell, R.A. Vicencio, M.I. Molina, and Yu.S. Kivshar, Phys. Rev. Lett. 97, 083901 (2006).
- (22) G. Kalosakas, Physica D 216, 44 (2006).
- (23) A.V. Savin, Y. Zolotaryuk and J.C. Eilbeck, Physica D 138, 267 (2000).
- (24) T. Ahn, R.S. MacKay, and J.-A. Sepulchre, Nonlinear Dynamics 25, 157 (2001); R.S. MacKay and J.-A. Sepulchre, J. Phys. A: Math. Gen. 35, 3985 (2002); J.-A. Sepulchre, in Localization and Energy Transfer in Nonlinear Systems, Proceedings of the Third Conference, San Lorenzo de El Escorial Madrid, edited by L. Vázquez, R.S. MacKay, and M.P. Zorzano (World Scientific, Singapore, 2003), p.102; M. Kastner and J.-A. Sepulchre, Disc. Cont. Dyn. Syst. B 5, 719 (2005).