Nonlinear surface modes and Tamm states in periodic photonic structures

Available online at www.sciencedirect.com

Wave Motion 45 (2007) 59–67 www.elsevier.com/locate/wavemoti

Nonlinear surface modes and Tamm states in periodic photonic structures Yuri S. Kivshar a

a,*

, Mario I. Molina

b

Nonlinear Physics Center, Research School of Physical Sciences and Engineering, Australian National University, Canberra, ACT 0200, Australia b Departmento de Fı´sica, Facultad de Ciencias, Universidad de Chile, Casilla 653, Santiago, Chile Received 15 November 2006; received in revised form 11 April 2007; accepted 12 April 2007 Available online 29 April 2007

Abstract We study the formation of nonlinear localized modes and discrete surface solitons near the edges or interfaces of weakly coupled nonlinear optical waveguides, one-dimensional photonic crystals. We draw an analogy between the staggered nonlinear surface optical modes and the surface Tamm states known in the electronic theory. We discuss the crossover between discrete solitons inside the array and surface solitons at the edge of the array by analyzing the families of even and odd nonlinear localized modes located at finite distances from the edge of a waveguide array. Then, we study the formation of guided modes localized at an interface separating two different periodic photonic lattices. Employing the effective discrete model, we analyze linear and nonlinear interface modes and also predict the existence of stable interface solitons including the hybrid staggered/unstaggered lattice solitons with the tails belonging to the spectral gaps of different types. Finally, we discuss briefly the recent experimental observation of discrete surface solitons and nonlinear Tamm states.  2007 Elsevier B.V. All rights reserved. Keywords: Surface waves; Tamm states; Surface solitons; Discrete solitons

1. Introduction Surface modes are a special type of waves localized at an interface between two different media. Surface states have been studied in different fields of physics, including optics [1,2] and nonlinear lattices [3]. In periodic systems, staggered modes localized at surfaces are known as Tamm states [4], first found as localized electronic states at the edge of a truncated periodic potential. Recently, it was predicted theoretically and demonstrated experimentally that nonlinear self-trapping of light near the edge of a waveguide array with self-focusing nonlinearity can lead to the formation of discrete surface solitons [5,6]. It was found that the self-trapped surface modes acquire some novel properties different from those of the discrete solitons in infinite lattices: discrete surface states can only exist above a certain *

Corresponding author. E-mail address: [email protected] (Y.S. Kivshar).

0165-2125/$ - see front matter  2007 Elsevier B.V. All rights reserved. doi:10.1016/j.wavemoti.2007.04.008

60

Y.S. Kivshar, M.I. Molina / Wave Motion 45 (2007) 59–67

threshold power and, for the same value of the power, up to two different surface modes can exist simultaneously. This can be understood as discrete optical solitons [7] localized near the surface but experiencing a repulsive force from the surface [8]. One of the important generalizations of these ideas is the concept of multi-gap surface solitons, i.e. mutually trapped surface states with the components associated with different spectral gaps [9]. Also, similar effects of light localization have been predicted and observed for self-defocusing nonlinear media together with the formation of surface gap solitons [10,11]. In this Letter, we overview our recent results on the physics of nonlinear surface localized modes and surface solitons. We consider a simple discrete nonlinear model described by the discrete nonlinear Schro¨dinger (NLS) equation and discuss two distinct cases: (a) surface solitons in a semi-infinite waveguide array and (b) interface solitons localized at the boundary between two different periodic photonic lattices. In particular, we explain the physical mechanism of the nonlinearity-induced stabilization of discrete surface modes and their existence above a certain power threshold. We also study the important generalization of the concept of nonlinear surface modes by analyzing nonlinear guided modes localized at an interface separating two different semi-infinite periodic photonic lattices. In the framework of an effective discrete model, we demonstrate that the analysis of linear interface states in such composite arrays provides an important tool for analyzing the interface solitons and their basic properties. We then find numerically the families of stable interface lattice solitons including a novel class of hybrid staggered/unstaggered lattice solitons with the tails localized in the spectral gaps of different types. We end with a brief discussion of recent experimental observation of discrete surface solitons and nonlinear Tamm states. 2. Surface solitons in waveguide arrays We study a semi-infinite array of identical, weakly coupled nonlinear optical waveguides [as shown in the inset of Fig. 1a] described by the system of coupled-mode equations [13,7] for the normalized mode amplitudes En dE1 2 þ E2 þ cjE1 j E1 ¼ 0; dz dEn 2 þ ðEnþ1 þ En1 Þ þ c jEn j En ¼ 0; i dz

i

ð1Þ

where n P 2, the propagation coordinate z is normalized to the intersite coupling V, En are defined in terms of 1=2 the actual electric fields En as En ¼ ð2V k0 g0 =pn0 n2 Þ En , where k0 is the free-space wavelength, g0 is the freespace impedance, a is the normalized linear propagation constant of each waveguide, n2 and n0 are nonlinear and linear refractive indices of each waveguide, and c = ±1 defines focusing or defocusing nonlinearity, respectively.

Fig. 1. Examples of surface localized modes at b = 3 (left) and b = 3 (right) in the array of focusing (c = +1, left) and defocusing (c = 1, right) waveguides centered at different distances d = 0 (a), 1 (b), 2 (c), 3 (d) from the edge.

Y.S. Kivshar, M.I. Molina / Wave Motion 45 (2007) 59–67

61

We look for stationary modes of the waveguide array in the form En(z) = exp(ibz)En, where b is the nonlinearity-induced shift of the propagation constant. For c = 0, we use the ansatz En  sin(nk) and obtain the linear spectrum b = a + 2 cos k, (0 6 k 6 p), and no localized surface modes. The presence of nonlinearity in the model (1) can give rise to new localized states. To find those modes, we analyze the stationary Eq. (1) where, without loss of generality, we scale out the parameter a. For given b, the system of stationary equations is solved numerically by a multi-dimensional Newton– Raphson scheme. Since we are interested in surface localized modes, we look for the states with maxima near the surface that decay quickly away from the array edge. Similar to an infinite array, these states could be centered at a waveguide site, or centered between waveguides. In an infinite discrete chain, such modes are known as odd and even states, respectively. In our calculations, we take N = 51 waveguides and explore both focusing and defocusing nonlinearities looking for localized modes below and above the linear spectrum band, jbj < 2. Fig. 1a–d (left and right) shows examples of the nonlinear localized states centered at different sites near the surface, for both focusing (c = +1, b = 3) and defocusing (c = 1, b = 3) nonlinearities, respectively. The surface state centered at the site n = 1 and shown in Fig. 1(a) was predicted earlier by Makris et al. [5]. The existence of multiple localized states near the surface and their stability are important characteristics of an interplay between nonlinearity and discreteness of the array, on one hand, and the surface created by the lattice truncation, on the other. In both the cases, the states (b,c) describe a crossover regime between the modes (a) with the maximum amplitude at the surface and the modes (d) which are weakly affected by the presence of the surface. To analyze the linear stability of each nonlinear stationary state found numerically, we introduce a weak perturbation as En(z) = En + [un(z) + ivn(z)] exp(ibz), and obtain linear evolution equations for un and vn, that can be expressed in a compact form by defining the real vectors dU{un} and dV = {vn}, and real matrices A = {Anm} and B = {Bnm}, 2

Anm ¼ dn;mþ1 þ dn;m1 þ ðb þ 3cjEn j Þ dn;m ; 2

Bnm ¼ dn;mþ1 þ dn;m1 þ ðb þ cjEn j Þ dn;m :

ð2Þ

€ þ BA dU ¼ 0, With the above definitions, the combined linear equations can be written in the form, dU € þ AB dV ¼ 0, where the dot stands for the derivative in z. Therefore, linear stability of nonlinear localized dV modes is defined by the eigenvalue spectra of the matrices AB and BA. If any of the real eigenvalues is negative, the corresponding nonlinear stationary solution is unstable; otherwise, the solution is stable. Results of this analysis are consistent with the so-called Vakhitov–Kolokolov stability criterion P of nonlinear localized modes, and the solitons determined by the slope of the power dependence P ¼ n jEn j2 , i.e. the states with dP/db < 0 for b > 0 or dP/db > 0 for b < 0, should be unstable. Fig. 2 shows the power P of the localized surface states vs. the propagation constant for the modes in the focusing waveguides shown in Fig. 1a–d, and the corresponding curves for the modes of the defocusing waveguides are mirror images. Direct numerical simulations and stability analysis confirm the validity of the Vakhitov–Kolokolov stability criterion; the instability region decreases as the center of the localized mode gets shifted away from the array edge. Similarly, we have also found even localized modes, akin to the modes found earlier for a semi-infinite nonlinear lattice [3], and verified that all in-phase even modes, for the focusing nonlinearity, and out-of-phase odd modes, for defocusing nonlinearity are all unstable, similar to the case of an infinite array. In order to get a deeper insight into the physics of the surface modes, we calP of the nonlinear stabilization P 4 culate the effective energy of the mode H ¼  Pn ðEn Enþ1 þ En Enþ1 Þ  12 n jEn j as a function of the distance of 2 the collective coordinate of the mode X ¼ P 1 n njEn j from the surface, similar to the case of a defect [14]. We apply a constraint method and start from the solution centered at the site n for given values of b and P. Our goal is to obtain all intermediate solutions between the odd and even stationary configurations for the same power. We proceed as follows: (i) we calculate an odd stationary mode centered at n and obtain all {En} and the power P, (ii) fix the amplitude at the site n + 1 to be En+1 + , (iii) solve the Newton–Raphson equations for all remaining Em (m 5 n + 1) with the constraint that the power be kept at P, arriving at an intermediate state centered between n and n + 1, and finally (iv) vary  and repeat the procedure until we reach the even configuration, where the amplitudes at the sites n and n + 1 coincide.

62

Y.S. Kivshar, M.I. Molina / Wave Motion 45 (2007) 59–67

(a)

–4.5 –5 –5.5 –6

1 1.5 2 2.5 3 3.5 4 4.5

X

EFFECTIVE ENERGY

EFFECTIVE ENERGY

Fig. 2. Normalized power vs. propagation constant b for the surface modes shown in Fig. 1 located at different distances d = 0, 1, 2, 3 from the surface. Black curve corresponds to the discrete soliton in an infinite array.

–8.5

(b)

–9 –9.5 –10 1 1.5 2 2.5 3 3.5 4 4.5

X

Fig. 3. Effective energy of surface modes vs. coordinate X near the edge of the array: (a) below (P = 2.85) and (b) above (P = 4.05) threshold. Black dots correspond to the stationary solutions found without constraint.

Fig. 3a and b shows the effective energy of the surface mode in a semi-infinite array, Ueff(X)  H(X), calculated for two different power values. The extremal points of this curve defined by the condition dH/dX = 0 correspond to the stationary localized solutions in the system. In comparison with an infinite array, the truncation of the waveguide array introduces an effective repulsive potential, that is combined with the periodic (Peierls–Nabarro) potential of an infinite waveguide array. As a result, discrete surface modes are possible neither in the linear regime nor in the continuous limit. As we see from Fig. 3a, for low powers there exists no solution of the equation dH/dX = 0 at the surface site n = 1; this corresponds to the fact that no surface state is found below the power threshold [5]. However, the modes localized at the sites n P 2 are still possible. If the power exceeds the threshold P = 3.26, discreteness overcomes a repulsive force of the surface and the surface localized state becomes possible, as shown in Fig. 3b. The correspondence between the stationary solutions found without constraints (black dots) and the solutions obtained as extremal points using the constraint methods is perfect. As expected, all odd modes are stable compared to even modes, and they all correspond to the condition dH/dX = 0. 3. Interface localized modes and hybrid lattice solitons In this section, we study another important concept in the theory of nonlinear surface modes. We analyze linear and nonlinear optical guided modes localized at an interface separating two different semi-infinite periodic photonic lattices. In the framework of an effective discrete model we demonstrate that the analysis of lin-

Y.S. Kivshar, M.I. Molina / Wave Motion 45 (2007) 59–67

63

ear interface states in such composite arrays provides an important tool for analyzing the interface solitons and their basic properties. We then find numerically the families of stable interface lattice solitons including a novel class of hybrid staggered/unstaggered lattice solitons with the tails localized in the spectral gaps of different types. We consider an interface separating two different semi-infinite arrays of optical waveguides (as shown in the top of Fig. 4) described by the system of coupled-mode equations [7] for the normalized mode amplitudes En: i

dEn 2 þ n En þ ðEnþ1 þ En1 Þ þ cjEn j En ¼ 0; dz

ð3Þ

where the waveguide interface is introduced by the condition: 0 at n = 0, and n = A or n = B for negative or positive n, respectively. jnj First, we look for linear surface modes in the form En ¼ An expðibzÞ localized near the interface waveguide with n = 0, and obtain ffi the condition n+/n = A0/B0 and the dispersion relation b ¼ 0 þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ð þ  Þð1  1 þ 4=A0 B0 Þ, where A0  A  0 and B0  B  0. Fig. 4 summarizes our results for A0 B0 2 the existence of such localized states on the parameter plane (A0, B0), as well as displays examples of localized modes corresponding to different existence regions. We note the existence of two sectors where no localized states exist (shaded regions). One of them is bounded by the curves B0 = 4/jA0j (for 1 < A0 < 0) and B0 = A0/(1  A0) (for 1 < A0 < 1). The other one is obtained by a reflection across the origin. Inside these regions either jn+1j, jn1j or both exceed one. One of the important observations that follows from our analysis is the existence of hybrid staggered/ unstaggered interface modes for the opposite sign of the propagation constant mismatches of two lattices. These modes have the tails localized in the bandgaps of different types, i.e. above (for one array) and below (for the other array) of the spectral band. In Fig. 5a and b we show the propagation constant b of the localized modes as a function of the interface parameter 0, for the characteristic cases (a) A = 0.6, B = 0.6, and (b) A = 3, B = 3. We note that the mode always lies outside the linear spectral bands, but its structure depends strongly on the overlap of the

B0 4

33 22 11 –44

–22

22

4 A0

–1 –2 –3 –4

Fig. 4. Phase diagram of different types of localized interface modes. No localized modes exist inside the shaded areas. The insets show examples of localized modes corresponding to different values of A0  A  0 and B0  B  0. Top: schematic structure of the waveguide interface.

Y.S. Kivshar, M.I. Molina / Wave Motion 45 (2007) 59–67 10

10

5

5

0

0

–5

–5

β

64

(a) –4

–2

0

2

(b)

4

–4

–2

0

0

2

4

0

Fig. 5. Families of localized modes shown as the dependencies b vs. 0, for the cases: (a) A = 0.6, B = 0.6 and (b) A = 3, B = 3. No localized modes exist inside the shaded regions. The dashed lines mark the spectral bands of both arrays, which in the case (b) do not overlap with each other.

bands, so that the hybrid modes appear for a relatively large band mismatch, as shown in Fig. 2b (middle curve). The analysis of linear localized interface modes in such an array provides an important information about the existence of nonlinear interface modes – lattice surface solitons. Next, we consider two semi-infinite nonlinear waveguide arrays characterized by the propagation constants A and B that are joined by an interface waveguide with the propagation constant 0. We focus on the interface modes defined by having their centers at either the first of the A waveguides or the first of the B waveguides. We find different classes of nonlinear localized modes and study their linear stability numerically via a multidimensional Newton–Raphson method. First, we consider the case 0 = A = B = 0.6 and c = +1. In the linear limit, this corresponds to a negative vertical line in Fig. 1 where no localized modes exist. The presence of nonlinearity shifts the propagation constant of the mode to the left, until it gives rise to an unstaggered localized mode. Therefore, we predict the lowest nonlinear interface mode to be unstaggered. This is indeed confirmed by our numerical computations, and the family of the lowest-order interface nonlinear modes is shown in Fig. 6, where the upper/lower branch corresponds to unstable/stable modes. Next, we consider the case A = 3, 0 = 0 and B = 3. Results are summarized in Fig. 7 which shows the dependence of the power vs. propagation constant for several low-power modes. We note that the lowest mode extends all the way to zero power, and therefore it corresponds in that limit to the linear mode induced by three concurrent dissimilar propagation constants. More importantly, the mode amplitudes show now a hybrid character, being unstaggered in one side of the interface and staggered on the other. In addition, we find many other types of nonlinear interface modes including twisted and flat-top modes, as well as the modes with different location of their centers relative to the interface, as discussed earlier for a semiinfinite waveguide array [8]. As a special limit of those modes, we find also the kink surface modes which are extended in one direction being localized in the other. 1.5 4.5

b power

a

1 0.5

4 –4 –2 1.5

3.5

0

2

4

b

1

a

6

0.5

3 3.4

3.6

3.8

4

–4 –2

0

2

4

6

waveguide

Fig. 6. Power vs. propagation constant for the interface unstaggered solitons centered on the first A waveguide (for A = 0 = 0.6, B = 0.6). The solid (dashed) curve refers to the stable (unstable) branches. Insets show two examples of the interface modes.

Y.S. Kivshar, M.I. Molina / Wave Motion 45 (2007) 59–67

a

1

20

65

0

power

15 –1 10

–4 –2

0

2

4

2

b a

5

b

0 –2

0

0.5

1

1.5

–4 –2 0

2

4

waveguide Fig. 7. Left: power vs. propagation constant for the hybrid interface staggered/unstaggered lattice solitons for A = 3, 0 = 0, B = 3. Solid (dashed) curves refer to the stable (unstable) branches. Right: examples of two lowest-lying hybrid interface solitons.

10

a

0.6

10

8

8

6

6

4

4

2

2

b

z

efficiency

0.8

a

0.4 0.2

b

0 1

2

3

4

5

6

7

8

–4 –2

0

2

4

–4 –2

0

2

4

input power Fig. 8. Left: trapping efficiency for the generation of interface solitons. The red (blue) curve denotes the case when the initial input is centered on the first A (B) waveguide. Right: evolution of initial states marked ‘a’ and ‘b’ (for A = 0.6, B = 0.6). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Next, we study the evolution of all types of nonlinear interface modes checking in this way our stability results and analyzing the evolution scenario for the unstable modes. In all cases examined, we observe that the unstable modes decay into the unstaggered fundamental mode by emitting radiation. In particular, for the unstaggered unstable modes, this decay is rapid and it proceeds vertically. For other cases, the unstable localized modes decay into higher-power generalizations of the fundamental mode. Also, we observe that unstable unstaggered modes decay much faster than unstable twisted, flat-top, or dark-like nonlinear modes. Finally, we study how an interface lattice soliton can be generated in experiment by exciting a single waveguide which is either the very first of the A waveguides, or the very first of the B waveguides. We define the trapping efficiency as the power fraction Pout/Pin remaining in the 10 central waveguides. Results for the trapping efficiency are shown in Fig. 5(top), where we have used 201 waveguides, a total evolution length of 20. The bottom portion of Fig. 8 shows the evolution of the initial states marked by the points ‘a’ and ‘b’ in Fig. 5(top). We would like to emphasize that the results obtained here can be easily generalized to the case of surface gap solitons which were predicted theoretically [10] and observed experimentally [11] in periodic photonic lattices with defocusing nonlinearity, where surface solitons appear in the gaps of the photonic bandgap spectra or their overlaps. 4. Experimental observations of surface optical solitons Shortly after the first prediction of formation of discrete optical solitons near the boundary of a self-focusing cubic nonlinear waveguide array [5], the first observation of these surface solitons was carried out [6]. The system consisted of a 101 AlGaAs waveguide array, into which a 1 kHz train of 1 ps pulses at 1550 nm from an optical amplifier, was focused onto the entrance facet of the waveguide array and at near the first channel. The input intensity distribution was shaped by a lens to accommodate the profile of the fundamental mode of

66

Y.S. Kivshar, M.I. Molina / Wave Motion 45 (2007) 59–67

Fig. 9. Theoretical prediction (a, b) and experimental observation (c, d) of nonlinear Tamm states in a truncated photonic lattice. (a) Schematic of the waveguide array geometry; (b) theoretical profile of a nonlinear Tamm statea surface gap soliton [4]; (c) threedimensional representation of the nonlinear surface state observed above the localization threshold [3]; (x, y) are the horizontal and vertical sample coordinates, respectively. (d) Experimental plane–wave interferogram demonstrating the staggered phase structure of the nonlinear Tamm state (from Ref. [11]).

a single channel. At low beam powers, discrete diffraction was observed, while at higher powers, a partial collapse of the diffraction was observed, until at high enough power the formation of a discrete surface soliton was observed [6]. While these surface optical solitons can be understood in terms of selfrapping of light at large enough powers, the first observations of bona fide nonlinear Tamm states came with the observation of surface gap solitons at the edge of a self-defocusing LiNbO3 array [11]. In this experiment, the array consisted of 100 single-mode optical waveguides fabricated by titanium in-difussion process in a monocrystal x-cut lithium Niobate wafer. The LiNbO3 exhibits a strong photovoltaic effect which leads to defocusing saturable nonlinearity at visible wavelengths. An extraordinarily polarized probe beam from a cw Nd:YVO4 laser (k = 532 nm) is focused by a microscope objective (20·) to a full width at half maximum (FWHM) of 2.7 lm at the input face of the sample, and injected into the waveguide at the edge of the array. The propagated wave packet at the output of the sample was imaged onto a CCD camera. The FWHM of the individual waveguide mode is 6 lm and 3 lm in horizontal and vertical directions, respectively, allowing for a single-waveguide input coupling. The waveguide array was externally illuminated by a white-light source in order to control the nonlinear response time. At low laser power (0.1 lW), two major effects were observed. First, due to coupling between neighboring waveguides the probe beam experiences discrete diffraction and spreads out in the horizontal plane upon propagation. Second, the beam shifts dramatically to the right indicating a strong repulsive effect of the surface. Increasing the laser power leads to spatial beam self-action through the defocusing photovoltaic nonlinearity. This nonlinearity-induced suppression of the surface repulsion leads to partial self-trapping at the surface, with a tail of intensity lobes extending into the periodic structure. A series of zero intensity points between these lobes indicates the self-induced dynamic formation of a staggered phase structure, or nonlinear Tamm state, as shown in Fig. 9. Similar experimental observation of staggered surface waves in a waveguide array in copper-doped lithium niobate substrate, was carried out in [12]. It must be noted here that, in all cases the experimental results and the theoretical predictions based on the effective discrete model were in good qualitative agreement. 5. Conclusions Employing the simple nonlinear model described by the discrete NLS equation, we have analyzed the physics and properties of different types of nonlinear localized modes that can be excited near the edge of a semiinfinite waveguide array. We have clarified the mechanism of the nonlinearity-induced stabilization and power threshold for the surface solitons. We have generalized this concept to the case of linear and nonlinear optical guided modes localized at the interface separating two different semi-infinite periodic photonic lattices. In the framework of an effective discrete nonlinear model, we have demonstrated the existence of stable interface lat-

Y.S. Kivshar, M.I. Molina / Wave Motion 45 (2007) 59–67

67

tice solitons including the hybrid staggered/unstaggered discrete solitons with the tails that belong to different spectral gaps. Our approach can be applied to other types of nonlinear discrete surface modes such as flat-top modes and twisted modes, as well as to the case of staggered modes in defocusing waveguides and two-dimensional surface solitons and their spatiotemporal generalizations. Acknowledgements We thank our colleagues C.R. Rosberg, D.N. Neshev, W. Krolikowski, A. Mitchell, and R.A. vicencio for useful collaboration in some of the problems presented in this paper. This work was supported by Conicyt and Fondecyt Grants 1050193 and 7050173 in Chile, and in part by the Australian Research Council. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]

P. Yeh, A. Yariv, A.Y. Cho, Optical surface waves in periodic layered media, Appl. Phys. Lett. 32 (1978) 104. W.J. Tomlinson, Surface wave at a nonlinear interface, Opt. Lett. 5 (1980) 323. Yu.S. Kivshar, F. Zhang, S. Takeno, Nonlinear surface modes in monoatomic and diatomic lattices, Physica D 113 (1998) 248. I.E. Tamm, A possible kind of electron binding on crystal surfaces, Z. Phys. 76 (1932) 849. K.G. Makris, S. Suntsov, D.N. Christodoulides, G.I. Stegeman, A. Hache, Discrete surface solitons, Opt. Lett. 30 (2005) 2466. S. Suntsov, K.G. Markis, D.N. Christodoulides, G.I. Stegeman, A. Hache´, R. Morandotti, H. Yang, G. Salamo, M. Sorel, Observation of Discrete Surface Solitons, Phys. Rev. Lett. 96 (2006) 063901. Yu.S. Kivshar, G.P. Agrawal, Optical solitons: from fibers to photonic crystals, Academic, San Diego, 2003. M. Molina, R. Vicencio, Yu.S. Kivshar, Discrete solitons and nonlinear surface modes in semi-infinite waveguide arrays, Opt. Lett. 31 (2006) 1693. I.L. Garanovich, A.A. Sukhorukov, Yu.S. Kivshar, M. Molina, Surface multi-gap vector solitons, Opt. Express 14 (2006) 4780. Ya.V. Kartashov, V.V. Vysloukh, L. Torner, Surface gap solitons, Phys. Rev. Lett. 96 (2006) 073901. C.R. Rosberg, D.N. Neshev, W. Krolikowski, A. Mitchell, R.A. vicencio, M.I. Molina, Yu.S. Kivshar, Observation of Surface Gap Solitons in Semi-Infinite Waveguide Arrays, Phys. Rev. Lett. 97 (2006) 083901. E. Smirnov, M. Stepic, C.E. Ruter, D. Kip, V. Shandarov, Observation of staggered surface solitary waves in one-dimensional waveguide arrays, Opt. Lett. 31 (2006) 2338. D.N. Christodoulides, R.I. Joseph, Discrete self-focusing in nonlinear arrays of coupled waveguides, Opt. Lett. 13 (1988) 794. Yu.S. Kivshar, F. Zhang, A.S. Kovalev, Stable nonlinear heavy-mass impurity modes, Phys. Rev. B55 (1997) 14265.