Aug 22, 2008 - arXiv:0806.1087v2 [nlin.PS] 22 Aug 2008 . Surface Solitons in Three Dimensions. Q.E. Hoq,1 R. Carretero-González,2 P.G. Kevrekidis,3 B.A..
.
Surface Solitons in Three Dimensions Q.E. Hoq,1 R. Carretero-Gonz´ alez,2 P.G. Kevrekidis,3 B.A. 4 5 Malomed, D.J. Frantzeskakis, Yu.V. Bludov,6 and V.V. Konotop6, 7
arXiv:0806.1087v2 [nlin.PS] 22 Aug 2008
1
Department of Mathematics, Western New England College, Springfield, MA, 01119, USA 2 Nonlinear Dynamical Systems Group, Department of Mathematics and Statistics, and Computational Science Research Center, San Diego State University, San Diego CA, 92182-7720, USA 3 Department of Mathematics and Statistics, University of Massachusetts, Amherst MA 01003-4515, USA 4 Department of Physical Electronics, Faculty of Engineering, Tel Aviv University, Tel Aviv 69978, Israel 5 Department of Physics, University of Athens, Panepistimiopolis, Zografos, Athens 157 84, Greece 6 Centro de F´ısica Te´ orica e Computacional, Universidade de Lisboa, Complexo Interdisciplinar, Avenida Professor Gama Pinto 2, Lisboa 1649-003, Portugal 7 Departamento de F´ısica, Faculdade de Ciˆencias, Universidade de Lisboa, Campo Grande, Ed. C8, Piso 6, Lisboa 1749-016, Portugal. (Dated: To appear in Phys. Rev. E, 2008) We study localized modes on the surface of a three-dimensional dynamical lattice. The stability of these structures on the surface is investigated and compared to that in the bulk of the lattice. Typically, the surface makes the stability region larger, an extreme example of that being the threesite “horseshoe”-shaped structure, which is always unstable in the bulk, while at the surface it is stable near the anti-continuum limit. We also examine effects of the surface on lattice vortices. For the vortex placed parallel to the surface this increased stability region feature is also observed, while the vortex cannot exist in a state normal to the surface. More sophisticated localized dynamical structures, such as five-site horseshoes and pyramids, are also considered.
I.
INTRODUCTION
Surface waves have been a subject of interest in a variety of contexts, including surface plasmons in conductors [1] and optical solitons in waveguide arrays [2] in physics, surface waves in isotropic magnetic gels [3] in chemistry, water waves in the ocean in geophysical hydrodynamics, and so on. Quite often, features exhibited by such wave modes have no analog in the corresponding bulk media, which makes their study especially relevant. In particular, a great deal of interest has been drawn to nonlinear surface waves in optics. It was shown theoretically [4] and observed experimentally [5] that discrete localized nonlinear waves can be supported at the edge of a semiinfinite array of nonlinear optical waveguide arrays. Such solitary waves were predicted to exist not only in selffocusing media (as in the above-mentioned works), but also between uniform and self-defocusing media [4, 6], or between self-focusing and self-defocusing media (e.g. in [7]). They have been subsequently observed in media with quadratic [8] and photorefractive [9, 10] nonlinearities. In the two-dimensional (2D) geometry, stable topological solitons have been predicted in a saturable medium [11], which constitute generalizations to lattice vortex solitons predicted in Ref. [12]. Quasi-discrete vortex solitons have been experimentally observed in a selffocusing bulk photorefractive medium [13]. Theoretical predictions for a variety of species of discrete 2D surface solitons [14, 15, 16, 17, 18], corner modes [15, 17], as well as surface breathers [17], were reported. Subsequent work has resulted in the experimental observation of 2D surface solitons, both fundamental ones and multi-
pulse states, in photorefractive media [19], as well as in asymmetric waveguide arrays written in fused silica [20]. Recently, surface solitons in more complex settings, such as chirped optical lattices in 1D and 2D [21, 22], at interfaces between photonic crystals and metamaterials [23], and in the case of nonlocal nonlinearity [24, 25], have emerged. Nearly all these efforts have been aimed at the study of surface solitons in 1D and 2D geometries. The only 3D setting examined thus far assumed a truncated bundle of fiber-like waveguides, incorporating the temporal dynamics in longitudinal direction to produce 3D “surface light bullets” in Ref. [26] (the respective 2D surface structures were examined in Ref. [27]). Our aim in the present work is to extend the analysis to surface solitons in genuine 3D lattices. Our setup is relevant to a variety of applications including, e.g., crystals built of microresonators trapping photons [28], polaritons [29], or Bose-Einstein condensates in the vicinity of an edge of a strong 3D optical lattice [30, 31]. In particular, we report results for discrete solitons at the surface of a 3D lattice, i.e., 3D localized states that are similar to relevant objects studied in the 2D setting of Ref. [14]. Thus, we will study localized states such as dipoles and “horseshoes” abutting on a set of three lattice sites, but also states that are specific to the 3D lattice. A variety of species of such solitons is examined below, and their stability on the surface is compared to that in the bulk. Some localized structures, such as dipoles, may be placed either normal or parallel to the surface. We demonstrate that, typically, the enhanced contact with the surface increases the stability region of the struc-
2 ture. Physically, this conclusion may be understood by the fact that the surface reduces the local interactions to fewer neighbors, rendering the system “more discrete”, hence more stable (by pushing the medium further away from the continuum limit, where all solitons would be unstable against the collapse). This effect is remarkable, e.g., for the three-site horseshoes which are never stable in the bulk, but get stabilized in the presence of the surface. However, the surface may also have an adverse effect, inhibiting the existence of a particular mode. The latter trend is exemplified by discrete vortices, which, if placed parallel to the surface, feature enhanced stability as compared to the bulk-mode counterpart, but cannot exist with the orientation perpendicular to the surface. Surface-induced effects of a different kind, which are less specific to discrete systems, are induced by the interaction of a particular localized mode with its fictitious “mirror image”. In terms of lattice models, the approach based on the analysis of the interaction of a real mode with its image was proposed in Ref. [32]. To formulate the model, we introduce unit vectors e1 = (1, 0, 0), e2 = (0, 1, 0), and e3 = (0, 0, 1) and deP3 fine lattice sites by n = j=1 nj ej with integer nj . We assume that the lattice occupies a semi-infinite space, n3 ≥ 1, and its dynamics obeys the discrete nonlinear Schr¨odinger (DNLS) equation in its usual form, iφ˙ n + ε∆φn + σ|φn |2 φn = 0.
(1)
Here φn is a complex discrete field, ε is the coupling constant, φ˙ n stands for the time derivative, the parameter σ = ±1 determines the sign of the nonlinearity (focusing or defocusing respectively), and ∆φn is the 3D discrete Laplacian: 3 X φn+ej + φn−ej − 2sφn , ∆φn≡
(2)
j=1
for n3 ≥ 2, while for n3 = 1 the term with subscript index n − e3 is to be dropped (note that e3 is the direction normal to the surface). It is interesting to point out here that an approach towards understanding the dynamics of Eq. (1) in the vicinity of the surface can be based on the above-mentioned concept of the fictitious mirror image, formally extends the range of n3 up to n3 = −∞, supplementing the equation with the anti-symmetry condition, φn1 ,n2 ,−n3 ≡ φn1 ,n2 n3 .
(3)
Indeed, this condition implies φn1 ,n2 ,0 ≡ 0, which is equivalent to the summation restriction in Eq. (2) as defined above. To confine the analysis to localized solitary wave modes, we impose zero boundary conditions, φn → 0 at n1,2 → ±∞ and n3 → ∞. Additionally, s = ±1 in Eq. (2) —this parameter is introduced for convenience (see Sec. III B) and can be freely rescaled using the transformation φ → φ eiνt for an appropriate choice of ν and
time rescaling. Stationary solutions to Eq. (1) will be sought for as φn = exp (iΛt)un , where Λ is the frequency and the lattice field un obeys the equation (Λ − σ|un |2 )un − ε∆un = 0.
(4)
Our presentation is structured as follows. The following section recapitulates the necessary background for the prediction of the existence and stability of lattice solitons. In section III, we report a bifurcation analysis for various surface states, treated as functions of coupling constant ε, with emphasis on the comparison with bulk counterparts of these states. Section IV reports the study of the evolution of unstable surface states. Finally, section V summarizes our findings and presents our conclusions. II.
THE THEORETICAL BACKGROUND
First, we outline some general properties of the model. Equation (1) conserves two dynamical invariants, namely the norm N , N=
∞ X
n3 =1 n1,2 =−∞
|φn |2 ,
(5)
and the Hamiltonian H, ∞ 3 X X σ ∗ ε φn (φn+ej − sφn ) + c.c. + |φn |4 , H= 2 n =1 j=1 3
n1,2 =−∞
(6) where the asterisk stands for complex conjugation. Stationary solutions to Eq. (4) with σ = ±1 are connected by the staggering transformation [17, 33]: if un is a solution for some Λ and σ = +1, then (−1)n1 +n2 +n3 un is a e = 12s − Λ and σ = −1. Consequently, it is solution for Λ sufficient to perform the analysis of stationary solutions, including their stability, for a single sign of the nonlinearity; thus, below we will fix σ = +1 (corresponding to the case of onsite self-attraction). Solutions to Eq. (4) in half-space n3 ≥ 1, subject to boundary condition φn = 0 for n3 = 0, as defined above, may be continued anti-symmetrically for the entire 3D space by setting Un ≡ un for n3 ≥ 1 and Un ≡ −un for n3 ≤ −1. Then, according to results of Ref. [34], this leads to an immediate conclusion, namely that there exists a minimum norm Nmin necessary for the existence of localized surface states in the present model. In other words, no surface modes survive in the limit of N → 0. In this connection, it is relevant to note that numerical findings that will be presented below were obtained, of course, for finite cubic lattices where, strictly speaking, there is no lower limit for N necessary for the existence of localized modes [17]. At this point, we have to specify that speaking about localized modes in a finite lattice we understand solutions which are localized on a number of cites much smaller than the total number of sites
3 in the chosen direction used for numerical simulations. Next we recall that generally speaking, there exist several branches of the nonlinear localized modes, i.e. for a given ε one can find localized excitations at different values of the norm N . Using the natural terminology we refer to higher/lower branches speaking about solutions with larger/smaller norm. In this classification the surface modes we are dealing with correspond to higher branches of the solutions of the respective finite lattices, i.e., their norm cannot be made arbitrarily small (see also the relevant discussion below in Section III B). To find solution families, we start with the anticontinuum (AC) limit, ε = 0 [35]. In this limit, the lattice field is assumed to take nonzero values only at a few (“excited”) sites, which determines the profile of the configuration to be seeded. The continuation of the structure to ε > 0 is determined by the Lyapunov’s reduction theorem [36]. More specifically, the solution is expanded as a power series in ε, the solvability condition at each order being that the respective projection to the kernel generated by the previous order does not give rise to secular terms [35]. The linear stability is then studied, starting from the usual form of the perturbed solution, ∗
φn = eiΛt (un + δan eλt + δbn eλ t ),
(7)
where δ is a formal small parameter, and λ is a stability eigenvalue associated with eigenvector ψ = {an , b∗n }. Substituting this into Eq. (1) yields the linearized system, iλan = −ε∆an + Λan − 2|un |2 an − u2n b∗n , −iλb∗n = −ε∆b∗n + Λb∗n − 2|un |2 b∗n − u∗n 2 an , which can be cast in the form ! ! H(1,1) H(1,2) A H(2,1)
H(2,2)
B
= iλ
A B
!
,
(8)
where A and B are vectors composed by elements an and b∗n , respectively, while the matrices H(p,q) (p, q ∈ {1, 2}) are given by, (1,1)
(2,2)
Hn,n′ = −Hn,n′ = δn,n′ Λ + 6sε − 2|un′ |2 −ε (1,2)
3 X j=1
δn+ej ,n′ + δn−ej ,n′ ,
(9)
(2,1) ∗
Hn,n′ = −Hn,n′ = −δn,n′ u2n′ .
An underlying stationary solution is (spectrally) unstable if there exists a solution to Eq. (8) with Re(λ) > 0. Otherwise, the stationary solution is classified as a spectrally stable one. As explained in Ref. [37], the Jacobian of the above mentioned solvability conditions is intimately connected to the full eigenvalue problem. More specifically, if the eigenvalues γ of the M × M eigenvalue problem of the Jacobian (where M is the number of excited sites at
the AC limit), then the near-zero eigenvalues of√the full stability problem can be predicted to be λ = 2γεp/2 , where p is the number of lattice sites that separate the adjacent excited nodes of the configuration at the AC limit.
III. A.
THE BIFURCATION ANALYSIS
Existence and stability of surface structures
In this section we study, by means of numerical methods, the existence and stability of various 3D configurations and compare the results to the corresponding analytical predictions. These configurations are obtained by starting from the AC limit (ε = 0), and are continued to ε > 0, using fixed-point iterations. For all the numerical results presented in this work, we fix the normalization Λ = 1 [see Eq. (4)], and use a lattice of size 13× 13× 13, unless stated otherwise. Also, for the presentation of the numerical results, we replace the triplet (n1 , n2 , n3 ) by (l, n, m), i.e., the surface corresponds to m = 1. We start by examining dipoles aligned parallel or normal to the surface. Panel (a) in Fig. 1 shows the norm of such states versus coupling constant ε, while panel (b) depicts the imaginary part of the stability eigenvalues for the bulk dipole, produced by the theory outlined in the previous section [dashed (black) lines], and by the numerical computations [solid (blue) lines]. It is worth mentioning that, for all the different configurations that we report in this manuscript, we display the imaginary part of the stability eigenvalue only for the bulk mode since the difference between the curves for the different variants (bulk, parallel or normal to the surface) is minimal. It should be noted however that the contact with the surface may produce higher order (smaller) eigenvalues that are not present in its bulk counterpart (results not shown here). The theoretical prediction for the sta√ bility eigenvalues is λ = ±2 εi, which, as expected, is the same as in an out-of-phase (twisted ) 1D mode analyzed in Ref. [37], since the structure is nearly onedimensional, along the line connecting the two excited sites. Panel (c) in Fig. 1 compares the largest instability growth rate as a function of ε for the bulk dipoles [dashdotted (green) line] and those oriented normally and parallel to the surface [dashed (red) and solid (blue) lines, respectively]. It is seen that the stability interval of the dipoles increases as its contact with the surface strengthens, in accordance with the arguments presented above. In the case of the bulk dipole, the instability occurs for values of the coupling constant in between ε0 = 0.114 and ε1 = 0.115. From now on, when reporting computed instability thresholds, we will use the lower bound for ε (e.g., ε0 in the above example) with the understanding that we always used the same resolution in ε. For the dipole set normally to the surface, we observe the onset of instability at ε = 0.117, while for the parallel-oriented one at ε = 0.120. In panels (e)–(h) of Fig. 1 we also depict
5 (a) 4 3 2 0
N 5 0.05
0.15
0.2
0 1
(b)
0.5 0 0
0.05
0.1 ε
0.15
0.1 (c)
0.05 0 0
Max(λ r )
0.04
0.05
0.1 ε
0.15
0.2
ε
(e)
l
λi
6
2 n 4 1m
3
−0.01
(g)
0
λr 0.01
2 (h)
7 6
0 −2
8
λi
l
7
8
0.122
2 (f)
8
6
0.12
0 −2
8
n 6
1m2
3
−5
0
λr
5 −3
x 10
FIG. 1: (Color Online) Results for the dipoles oriented parallel and normal to the surface. (a) Norm N versus the lattice coupling constant, ε. (b) Imaginary part of the linear stability eigenvalue: solid (blue) and dashed (black) lines correspond, respectively, to numerically found and analytically predicted forms. (c) Real part of the critical (in)stability eigenvalue: the dashed (red) and solid (blue) lines depict the normal- and parallel-oriented dipoles, respectively, while the dash-dotted (green) line corresponds to the bulk dipole. (d) (In)stability eigenvalue for the parallel surface dipole placed at distances from the surface starting from zero and up to five lattice periods away (curves right to left). (e)-(g) Configurations and (f)-(h) respective spectral stability planes just above the instability threshold. The level contours, corresponding to Re(ul,n,m ) = ±0.5 max {ul,n,m } are shown, respectively, in dark gray (blue) and gray (red). The instability thresholds for the dipoles oriented parallel and normally to the surface are, respectively, ε = 0.117 and ε = 0.120. For comparison, the threshold for the bulk dipole is ε = 0.114.
the shapes of the normal and parallel dipoles, just below the instability threshold, along with their corresponding spectral stability planes. The stabilizing effects exerted by the surface depend, in a great measure, on the distance of the configuration from the surface, namely, the further away the configuration from the surface, the lesser the effect is. This property is clearly seen in panel (d) of Fig. 1, where we
(b)
0.05
0 0
l
0.118
0.4
0.1
0.15 ε
0.2
0.25
0.3
(c)
8 0.116
0.3
0.2
(d)
0.114
0.2 ε
0.5
0.02 0 0.112
0.1
0 0
0.2
Max(λ r )
Max(λ r )
0.1 ε
0.1
0.2 ε
(d)
7 6
8
0.3 5
λi
λi
1
10 (a)
λi
N
4
6
2 n 4 1m
3
0.4
(e)
0
−5 −0.02
0
λr 0.02
FIG. 2: (Color Online) The stability of the three-site “horseshoe”. Panels are similar to those in Fig. 1. Panel (c) compares the critical stability eigenvalue, as a function of the lattice coupling, ε, for the surface and bulk horseshoes [solid (blue) and dashed-dotted (green) lines, respectively]. The bulk horseshoe is always unstable (due to a purely real, higher-order eigenvalue), while the corresponding surface configurations have a stability region (the corresponding eigenvalue becomes imaginary in this case). Panels (d)-(e) correspond to the surface horseshoe just above the stability threshold of ε = 0.239.
plot the (in)stability eigenvalue as a function of the coupling for several values of the separation of the parallel dipole from the surface. The curves, from right to left, depict the results for the dipole set at the distance of 0, 1, ..., 5 sites away from the surface (0 sites refers to the dipole sitting on the surface). As the panel demonstrates, the stability interval is reduced as the dipole is pulled away from the surface, converging towards a bulk dipole. Let us now consider the “horseshoe” configurations, for which the presence of the surface is critical to their stability. In Fig. 2 we depict the properties of a threesite horseshoe, which actually is a truncated version of a quadrupole, cf. the 2D situation [14]. As before, panel (a) in Fig. 2 shows the norm versus ε, while panels (b) and (c) compare the stability of the bulk horseshoe (the dash-dotted line) and ones built near the surface (the solid line). While the bulk horseshoes are always unstable, similar to their 2D counterparts [14], the ones placed near the surface are stable at small ε, destabilizing at ε = 0.239. Panels (d)-(e) in Fig. 2 show the configuration for the coupling just below the instability threshold, along with the respective spectral plane. The analytical √ expressions for stable eigenvalues are λ = 0, λ = ±2 3εi, λ = O(ε2 ), cf. the expressions obtained in Ref. [14] for the 2D horseshoes.
5
10 5 0
0.05
0.1
0.15 ε
0.2
0.25
1
λi
(b)
0.5 0 0
(a)
6 4 0
0.3
0.05
0.1
0.15 ε
0.2
0.25
0.04
0.2 (c)
0.02
0.04
ε
0.06
0.08
0.1
0.02
0.04
ε
0.06
0.08
0.1
(b)
0.5 0 0
0.3
Max(λ r )
(c)
0.02
0 0
0.05
0.1
0.15 ε
(d)
λi
8 l
7 6 10 n
5
3 4 1 2m
0.2
4 (e) 2 0 −2 −4 −0.02
0.25
0.3
0 0.065
0 λ 0.02 r
FIG. 3: (Color Online) The stability for the five-site horseshoe at the surface. Panels are identical to those in Fig. 2. In this case, the stability threshold is at ε = 0.211, while for the bulk 5-site horseshoe it is ε = 0.205. Panels (d)-(e) depict the configuration and the respective linear stability spectrum just above the critical point of ε = 0.211.
Figure 3 illustrates the same features as before but for the five-site horseshoe. Unlike its three-site cousin, the bulk five-site horseshoe is stable up to a critical value of the coupling, ε = 0.205, while the surface variant has it stability region ε < 0.211. The eigenvalues of the linearization in this case can be computed similar to those for the three-site modes [14], as outlined above (cf. also Ref. [35]), which eventually yields λ = 3.8042εi, λ = 2.8284εi, λ = 2.3511εi, λ = O(ε2 ), and λ = 0, in good agreement with the corresponding numerical results, as shown in panel (b) Fig. 3. Next we consider the quadrupole configuration, see Fig. 4. The surface again exerts a stabilizing effect, albeit a weaker one, when the quadrupole is placed normally and parallel to the surface. In the bulk, the quadrupole loses stability at ε = 0.068, while the normal and parallel surface quadrupoles have stability thresholds, respectively, at ε = 0.070 and ε = 0.071. The analytical approximation for the stability eigenvalues√in this case are √ λ = 8εi (a double eigenvalue), λ = 2 εi, and a zero eigenvalue, which accurately capture the numerical findings depicted in panel (b) of Fig. 4. In Fig. 5 we present the results for four-site vortices. This configuration, in contrast to the previous ones, is described by a complex solution. In the AC limit, the vortex occupies the same excited sites as the above-mentioned quadrupole, but the phase profile, {0, π/2, π, 3π/2}, emulates that of the vortex of charge 1 [12, 35]. The bulk four-site vortex (which was discussed in Ref. [38]) loses its stability at ε = 0.438, while the
0.07
(d) 8 7 6 5 8 7 3 2 n6 5 1 m
0.075
ε
0.08
2 (e) λi
0.1
l
λi
1
Max(λ r )
8
(a)
N
N
15
0 −2 −2
λr
0
2 −3 x 10
FIG. 4: (Color Online) The stability of quadrupole modes. The layout is similar to that in Fig. 3. In panel (c), due to the close proximity of the thresholds, the close-up is shown for the critical stability eigenvalue versus the lattice coupling constant, ε, for the parallel and normal surface modes, and the bulk one [solid (blue) and dashed (red) lines, and the dash-dotted (green) line, respectively]. The threshold for the bulk mode is ε = 0.068, while for the normal and parallel quadrupoles it is, respectively, ε = 0.070 and ε = 0.071. As before, panels (d) and (e) show the configuration just above the instability threshold along with its corresponding spectralstability plane.
vortex parallel to the surface features an extended stability region, up to ε = 0.505. However, the surface in this case prohibits the existence of a vortex that would be oriented normally to the surface layer, similarly to what was found for 2D lattice vortices [14]. The simplest explanation for the complete absence of the solution normal to the surface, compared with that of an existing vortex waveform parallel to the surface can arguably be traced in the interaction of such vortices in the half-space with their fictitious image (if the domain is equivalently extended to the full space). In the case of the vortex parallel to the surface, the situation is tantamount to the vortex cube structures examined in [39, 40], for which it was established in [40] that the persistence conditions are satisfied (and, in fact, that such structures consisting of two out-of-phase vortices should be linearly stable close to the AC limit). On the other hand, for the case normal to the surface, by examining the relevant interactions it can be observed (at an appropriately high order) that the persistence conditions of [35, 37, 40] can not be satisfied and hence the structure can not be continued beyond the AC limit. That is why the structure can never be observed to exist irrespectively of the smallness of ε. The next species of stationary lattice solutions is a
6
20 (a)
(a)
N
N
10
0.6
λi 0.2
0.4 ε
0.6
0.6
10 λi
l
0.4 ε
6
n5
2 3 4 1 m
0.8
0.2
0.25
0.3
0.05
0.1
0.15 ε
0.2
0.25
0.3
0.05
0.1
0.15 ε
0.2
0.25
0.3
(d)
0
−10 −0.01
0.15 ε
0.5 (c) 0 0
(e)
0.1
(b)
0 0
Max(λ r )
0.2
(d)
0.05
0.5
0.8
0.5 (c)
7 6 5 4
0
1
(b)
0 0
0
λ 0.01 r
FIG. 5: (Color Online) The stability of the four-site vortex in the grid of size 11× 11× 11. The dash-dotted and solid lines show the bulk vortex and the one parallel to the surface, respectively. The layout is similar to that of the above figures. Instability in the bulk occurs at ε = 0.438, and in the parallel surface vortex at ε = 0.505. The vortex cannot exist with the orientation normal to the surface. Panels (d) and (e) show the parallel surface vortex just above the instability threshold of ε = 0.485. As in the previous figures, the level contours corresponding to Re(ul,n,m ) = ±0.5 max {ul,n,m } are shown, respectively, in dark gray (blue) and gray (red), while the complementary level contours, defined as Im(ul,n,m ) = ±0.5 max {ul,n,m }, are shown by light gray (green) and very light gray (yellow) hues, respectively.
pyramid-shaped structure, with characteristics displayed in Fig. 6, whose base is a rhombus composed of four sites. The remaining out-of-plane vertex site must have phase 0 or π, since the phase values π/2 and 3π/2 at this site do not produce a solution. The full set of pyramids (bulk, normal, parallel —see panels (d)–(f) of Fig. 6) is completely unstable, as seen in panel (c) of Fig. 6, the surface producing no stabilizing effect on it. This strong instability actually arises at the lowest order in the √ analytical eigenvalue calculations, which yield λ = 2 5εi, √ λ = 2 2εi, λ = 2ε, λ = 0, and λ = O(ε2 ), once again in very good agreement with the full numerical results of Fig. 6.
B.
5
0.8
0.5 0 0
Max(λ r )
0.4 ε
l
λi
1
0.2
Small-amplitude modes in a finite lattice
Since our numerical investigation of the surface modes uses a finite lattice, which allow the existence of smallamplitude modes (ones with the zero threshold in terms of the norm — cf. discussion in Sec. II), here we briefly consider the modes in a finite lattice having the smallamplitude limit. Our aim is to show that these modes belong to lower branches, as compared with the “nor-
(e)
8 7 6 8 7 n 6
l
0
8 6 7m
(f)
8 7 6
l
10
8 7 n 6
8 7 6
3 4 1 2 m
8 7 n 6
1
4 3 2 m
FIG. 6: (Color Online) The instability of pyramid-shaped structures. This configuration abuts on the base in the form of a rhombus, and includes the out-of-plane site with zero phase. Three variants of this configuration are displayed in panels (d)–(f): bulk, normal and parallel to the surface, respectively. The stability of the three different variants of the pyramid is essentially identical, all three of them being unstable.
mal” surface modes considered above. To this end, we concentrate on the lattice of size M× M× M lattice, subject to the zero boundary conditions, which imply that discrete Laplacian (2) is modified at surfaces nj = 1 and nj = M (j = 1, ..., 3) by setting the fields at sites n − ej and n + ej , respectively, equal to zero. For the sake of definiteness, we fix here s = −1 in Eq. (2). To determine the norm N of small-amplitude modes we follow Ref. [17], and look for a solution to Eq. (4) with the amplitude un and coupling constant ε being represented as series un = ǫu0,n + ǫ2 u2,n + O(ǫ3 ), 2
(10)
3
ε = ε0 + ǫ ε2 + O(ǫ ), p in powers of small parameter ǫ ≡ 8N/(M + 1)3 ≪ 1, which vanishes in the limit of the infinite lattice (M → ∞); in other words, small ǫ characterizes the “largeness” of the lattice. We focus on real solutions here. Substituting expansions (10) into Eq. (4) and gathering terms of the same order in ǫ, we rewrite Eq. (4) in the form of a set of equations: Λuj,n − ε0 ∆uj,n = Fj,n .
(11) 3
Here F0,n = 0, F2,n = Λ(ε2 /ε0 )u0,n + (u0,n ) , hence
7 Eq. (11) with j = 0 gives rise to a linear eigenmode, sin
j=1
πnj mj M +1
,
(12)
5 0 0
with the respective approximation for the lattice coupling constant, = Λ 6 + 2
3 X j=1
cos
πmj M +1
−1 ,
0.4
(13)
parameterized by vector m = (m1 , m2 , m3 ). At the same time, considering the solvability conditions for j = 2, which amounts to demanding the orthogonality of F2,n and u0,n , we obtain corrections to the coupling constants, (m)
ε = ε0
−
(m) 3 ǫ2 ε0 Y 3 + δmj ,(M+1)/2 . 64Λ j=1
(14)
It follows from Eq. (14) that each of the linear modes (12) is uniquely extended into a small-amplitude nonlinear one. These modes are characterized by the linear dependence of the norm on coupling constant ε: (m) 8Λ(M + 1)3 ε0 − ε (15) N (m) = (m) Q3 . ε0 3 + δ m ,(M+1)/2 j j=1
From Eq. (15) it follows that each mode, parameterized by vector m, exists when ε belongs to the inter(m) val 0 ≤ ε ≤ ε0 . The validity of approximation (15) is corroborated by the coincidence of analytical and numer(m) ical results in the vicinity of ε0 (as shown in Fig. 7), where these modes reaches their small-amplitude limit. Such a property of these modes differs considerably from the case of the surface modes which do not possess the small-amplitude limit and require some minimal value of the norm (for the normal dipole, depicted in Fig. 7 by dash-dotted line, the minimal norm is ≈ 1.262). Panel (b) in Fig. 7 shows that only the mode, parameterized (m) by vector m = (1, 1, 1), is stable for ε close to ε0 , while other modes are completely unstable.
IV.
DYNAMICS
In this section we examine the nonlinear evolution of the various configurations, displaying the results in a set of figures (see Figs. 8–12). In each case, the evolution is initiated at a value of the coupling ε taken beyond the instability threshold, and an initial small random perturbation is applied in order to expedite the onset of the instability. All the figures display the evolution of the instability at six different moments of time, starting at t = 0, and
max(λr)
(m)
ε0
(a)
m=(1,2,2)
10
N
(m)
u0,n =
3 Y
15
0.2 0 0
m=(1,1,1) m=(1,1,2) 0.02
0.04
ε
0.06
0.08 (b)
m=(1,2,2) m=(1,1,2)
0.02
m=(1,1,1) 0.04
ε
0.06
0.08
FIG. 7: (Color Online) Low-amplitude modes in a finite grid of size 9×9×9 with Λ = 1.0. (a) Norm N versus coupling constant ε for several modes whose low-amplitude limit is parameterized by vectors m, as calculated numerically and predicted by approximation (15) (solid and dashed lines, respectively). For comparison, dash-dotted line depicts the norm for surface normal dipole. (b) Real part of the critical (in)stability eigenvalue, calculated numerically.
ending at a time well beyond the point at which the instability manifests itself. All configurations that were predicted above to be unstable through nonzero real parts of the (in)stability eigenvalue λ indeed exhibit instability dynamics, which eventually results in a transition to a different configuration. In the case of the dipoles and horseshoes, Figs. 8–10 show a spontaneous transition to monopole patterns, i.e., ones centered around a single excited site. On the other hand, in the case of the vortices and pyramids shown in Figs. 11–12, a few sites may remain essentially excited at the end of the evolution. The monopole is, obviously, the most robust dynamical state in the lattice system, with the widest stability interval, in comparison with other discrete structures. This simplest state becomes unstable, for given Λ, only at values of the coupling constant ε ≈ Λ [38]. Another structure with a relatively wide stability region is the dipole (the more stable the wider the distance between its constituent sites [39]), consonant with the observation that some of the structures (especially ones with a large number of excited sites, such as vortices and pyramids) dynamically transform into dipoles. Generally speaking, the exact scenario of the nonlinear evolution and the finally established state depend on details of the initial perturbation. In the figures, each configuration is shown by iso-level contours of distinct hues. In particular, dark gray (blue) and gray (red) are iso-contours of the real part of the solutions, while the light gray (green) and very light gray (yellow) colors depict the imaginary part of the same solutions. A case that needs further consideration is that of the three-site horseshoe. As observed from the stability analysis presented in Fig. 2, this horseshoe in the bulk gives rise to a small unstable purely real eigenvalue for all val-
8
6
7
7
7
7
8
8
7 m
n6
7 6
1
2
6
7 m
7
8
n6
t=43
t=48
7
7
l
l
(b) t=0
7
7
8
7
3
6
1
2
7
3
6
7
7
n6
1
2 m
l
7
l
t=100
l
t=51
7
3
n6
1
2 m
1
7
3
7 m
n6
1
2
2 m
3 8 7 n 61
2
2 m
3
7
8
8
7 3 4 n 6 1 2m
4 6 1 2 3
7 3 4 n 6 1 2m
3
FIG. 8: (Color Online) The evolution of unstable dipoles: (a) a bulk dipole; (b) and (c) dipoles placed parallel and normally to the surface, respectively. In all the cases, the dipole is subject to oscillatory instability, which is responsible for the eventual concentration of most of the norm at a single site (i.e., the transition to a monopole). Parameters are Λ = 1, ε = 0.2, the lattice has a size of 13× 13× 13, and times are indicated in the panels. All iso-contour plots are defined as Re(ul,n,m ) = ±0.75 = Im(ul,n,m ), and the initial configurations were perturbed with random noise of amplitude 0.01. The coding for the iso-contours is as follows: dark gray (blue) and gray (red) colors pertain to iso-contours of the real part of the solutions, while the light gray (green) and very light gray (yellow) colors correspond to the iso-contours of the imaginary part.
FIG. 9: (Color Online) The evolution of the unstable threesite horseshoes: (a) bulk three-site horseshoe and (b) the horseshoe oriented normally to the surface. In both cases, the unstable horseshoe is subject to an oscillatory instability, which leads to the eventual concentration of most of the norm in a single-site structure. The iso-contours and parameters are the same as in Fig. 8 except that ε = 0.3. (a) t=0
8 7 6
l
3 8 7 n 61
7 3 4 n 6 1 2m
8 t=24
6
8
n6
6m
8
8
n6
2m
8
8
t=50
n6
6m
8
8
t=18
8 7 6
8
6
2
8
4
t=50
8 7 6 8
n6
2m
4
6
6
8
8 7 6
8 7 6 t=24
4
6
t=17
4
8 7 6
6
8 7 6
8
l
l
2
8 7 6
t=26
8 7 6 6
t=23
8 7 6 8
(b) t=0
8
t=20
8
8 7 6
t=23
ues of ε, see the lower green dashed-dotted curve in panel (c) of the figure. Despite the presence of this eigenvalue, the evolution of the unstable bulk three-site horseshoe is predominantly driven by the unstable complex eigenvalues, if any (in fact, for ε > 0.226, see the dashed-dotted (green) line of Fig. 2.(c)). A careful analysis of the instability corresponding to the small purely real eigenvalue for ε < 0.226 (i.e., before the complex eigenvalues become unstable) reveals that the corresponding dynamics amounts to an extremely weak exchange of the norm between the two in-phase excited sites (see Fig. 2). The
6
l
7
l
7
l
7
l
t=100
2 m
8
4 6 1 2 3 8 7 6
l
61
t=53
7 n 61
7
7
3 8 7
t=51
8
8
t=100
l
61
2
8 7 6
8 7 6
l
61
2
3 8 7
t=35
8 7 6
t=63
l
7
7
l
l
l
7 8
t=27
4 6 1 2 3
9 6 6 7 8
7 8 9 n 6 6 7m
8 7 6
t=48
t=39
7
8
7 8 9 n 6 6 7m
t=43
8 (c) t=0
8
7 8 9 n 6 6 7m
8 7 6 8
2 m
l
8 7 6
8
3
7
t=50
8
3
8
9 6 6 7 8
8 7 6
(b) t=0
t=49
7
6
7
t=29.5
l
6
8
9 6 6 7 8
8 7 6
l
n6
l
7
7
t=26
l
l
t=50
l
t=47
l
t=38
6
l
7
8
l
7
l
6
l
6
l
7
8
l
7
l
6
8 7 6
l
6
t=25
8 7 6
l
7
t=22
l
7
8 7 6
l
7
(a) t=0
l
t=30
l
7
t=29
l
l
(a) t=0
n6
6
6m
2
8
4
8 7 6 8
n6
2m
4
FIG. 10: (Color Online) The evolution of unstable five-site horseshoes: (a) the bulk horseshoe, and (b) the five-site horseshoe oriented normally to the surface. In both cases, the unstable horseshoe is subject to an oscillatory instability, which triggers the transition to a monopole. The iso-contours and parameters are the same as in Fig. 9.
9
(b) t=0
7
7 n6 6 m
t=42
l 6
7
6
t=44
1
2
6 7
3
6
t=46
2
7
3
n6
2 1 m
3
1
2
3
l 6
7
6 7
n6
2 1 m
3
7
n6
2 1 m
3
FIG. 11: (Color Online) The evolution of unstable vortices: (a) the bulk vortex for ε = 0.3 and (b) the vortex parallel to the surface, for ε = 0.6 and Λ = 1. The iso-contour plots are defined by Re(ul,n,m ) = ±1 = Im(ul,n,m ).
norm exchange is driven by the corresponding unstable eigenfunction, which looks like a dipole positioned at the two aforementioned in-phase sites. The difficulty in observing this evolution mode is explained by the fact that, in the course of the norm exchange, only ∼ 0.1% of the total norm is actually transferred between the two sites. Furthermore, as mentioned earlier, the corresponding small real eigenvalue is completely suppressed by the surface (see panel (c) in Fig. 2). It is worth noting that such stable three-site horseshoe surface structures may also be generated by the evolution of more complex unstable waveforms, such as the five-site pyramids placed normally to the surface, see the bottom panel in Fig. 12.
V.
7 6 5 87 6 t=82
7
l 6
6
t=100
7
l
7
1
7 8 6 m
(b) t=0
7
l 6
8
t=43
7
l
7
7 n6 6 m
CONCLUSIONS
In this work, we have investigated localized modes in the vicinity of a two-dimensional surface, in the framework of the three-dimensional DNLS equation, which is a prototypical model of nonlinear dynamical lattices. We have found that the surface may readily stabilize localized structures that are unstable in the bulk (such as three-site horseshoes), and, on the other hand, it may inhibit the formation of some other structures that exist in the bulk (such as vortices which are oriented normally to the surface, although ones parallel to the surface do exist
7 6 5 87 n6
(c) t=0 8 7 6 87
6 t=212 8 7 6 87 n6
l l
l
l 7
8
7 6 5 87 n6
8 6 7
7 8 6 m
3 1 2
t=20
t=66
7 6 5 87 6 t=182
7 6 5 87 6 t=400
3 1 2 m
3 1 2
3 1 2 m
7 6 5 87 n6
7 6 5 87 n6
t=30 8 7 6 87
6 t=222 8 7 6 87 n6
8 6 7
7 8 6 m
l
8
6
7 6 5 87 n6
8 6 7
l
7 n6 6 m
8
l 6
7
7
7
l 6
6 6
t=100
7
l
7
7
8
l
t=47
7
l
6 6
7 6 5 87 6 t=268
3 1 2
3 1 2 m
l
7
8
l
7
l
t=31
6 6
l
7
6
7 6 5 87 6 t=188
3 1 2
l
6
l
l
l 6
7 6 5 87 6 t=156
t=96
l
7
t=16
7 6 5 87 n6
t=156 8 7 6 87 6 t=400 8 7 6 87 n6
3 1 2
3 1 2 m
l
7
7
(a) t=0
3 1 2
l
t=27
l
t=24
l
(a) t=0
3 1 2 m
3 1 2 m
FIG. 12: (Color Online) The evolution of unstable pyramids. Panels (a), (b), and (c) display, respectively, the transformation of a bulk pyramid, and of ones oriented normally and parallel to the surface, for ε = 0.2.
and have their stability region; a qualitative explanation to these features was proposed, based on the analysis of the interaction of the vortical state with its “mirror image”). The most typical surface-induced effect is the expansion of the stability intervals for various solutions that exist in the bulk and survive in the presence of the surface. This feature may be attributed to the decrease, near the surface, of the number of neighbors to which excited sites couple, since the approach to the continuum limit, i.e., the strengthening of the linear couplings to the nearest neighbors, is responsible for the onset of the instability or disappearance of all the localized stationary states in the three-dimensional dynamical lattice. On the other hand, while the techniques elaborated in Refs. [35, 37, 40] for the analysis of localized states in
10 bulk lattices are quite useful in understanding the dominant stability properties of the solutions, the surface gives rise to specific effects, such as the stabilization of higherorder solutions or the suppression of some types of vortex
structures, which cannot be explained by these methods. Therefore, it would be very relevant to modify these techniques, which are based on the Lyapunov-Schmidt reductions, so as to take the presence of the surface into regard.
[1] W. L. Barnes, A. Dereux, and T. W. Ebbesen, Nature 424, 824-830 (2003). [2] D. N. Christodoulides, F. Lederer, and Y. Silberberg, Nature 424, 817-823 (2003). [3] S. Bohlius, H. R. Brand and H. Pleiner, Z. Phys. Chem. 220, 97-104 (2006). [4] K. G. Makris, S. Suntsov and D. N. Christodoulides, G. I. Stegeman, Opt. Lett. 30, 2466 (2005); M. I. Molina, R. A. Vicencio, and Yu. S. Kivshar, Opt. Lett. 31, 1693 (2006). [5] S. Suntsov, K. G. Makris, D. N. Christodoulides, G. I. Stegeman, A. Hache, R. Morandotti, H. Yang, G. Salamo, and M. Sorel, Phys. Rev. Lett. 96, 063901 (2006). [6] Y. V. Kartashov, V. A. Vysloukh, and L. Torner, Phys. Rev. Lett. 96, 073901 (2006). [7] D. L. Machacek, E. A. Foreman, Q. E. Hoq, P. G. Kevrekidis, A. Saxena, D. J. Frantzeskakis, and A. R. Bishop Phys. Rev. E 74, 036602 (2006). [8] G. Siviloglou, K. G. Makris, R. Iwanow, R. Schiek, D. N. Christodoulides, G. I. Stegeman, Y. Min, and W. Sohler, Opt. Exp. 14, 5508 (2006). [9] C. R. Rosberg, D. N. Neshev, W. Kr´ olikowski, A. Mitchell, R. A. Vicencio, M. I. Molina, and Yu. S. Kivshar, Phys. Rev. Lett. 97, 083901 (2006). [10] E. Smirnov, M. Stepi´c, C. E. R¨ uter, D. Kip, and V. Shandarov, Opt. Lett. 31, 2338 (2006). [11] Y. V. Kartashov, A. A. Egorov, V. A. Vysloukh, and L. Torner, Opt. Exp. 14, 4049 (2006). [12] B. A. Malomed and P. G. Kevrekidis, Phys. Rev. E. 64, 026601 (2001). [13] D. N. Neshev, T. J. Alexander, E. A. Ostrovskaya, Yu. S. Kivshar, H. Martin, I. Makasyuk, and Z. Chen, Opt. Phys. Rev. Lett. 92, 123903 (2004). [14] H. Susanto, P. G. Kevrekidis, B. A. Malomed, R. Carretero-Gonz´ alez, and D. J. Frantzeskakis, Phys. Rev. E. 75, 056605 (2007). [15] K. G. Makris, J. Hudock, D. N. Christodoulides, G. I. Stegeman, O. Manela, and M. Segev, Opt. Lett. 31, 2274 (2006). [16] Y. V. Kartashov and L. Torner, Opt. Lett. 31, 2172 (2006). [17] Yu. V. Bludov and V. V. Konotop, Phys. Rev. E 76, 046604 (2007). [18] R. A. Vicencio, S. Flach, M. I. Molina and Yu. S. Kivshar, Phys. Lett. A 364, 274 (2007). [19] X. Wang, A. Bezryadina, Z. Chen, K. G. Makris, D. N. Christodoulides, and G. Stegeman, Phys. Rev. Lett. 98, 123903 (2007). [20] A. Szameit, Y. V. Kartashov, F. Dreisow, T. Pertsch, S. Nolte, A. T¨ unnermann, and L. Torner, Phys. Rev. Lett.
98, 173903 (2007). [21] Y. V. Kartashov, V. A. Vysloukh, and L. Torner, Phys. Rev. A 76, 013831 (2007). [22] M. I. Molina, Y. V. Kartashov, L. Torner, Yu. S. Kivshar, arXiv:0712.3179. [23] A. Namdar, I. V. Shadrivov, and Yu. S. Kivshar, Phys. Rev. A 75, 053812 (2007). [24] B. Alfassi, C. Rotschild, O. Manela, M. Segev, and D. N. Christodoulides, Phys. Rev. Lett. 98, 213901 (2007). [25] F. Ye, Y. V. Kartashov, and L. Torner, arXiv:0802.2521. [26] D. Mihalache, D. Mazilu, F. Lederer and Yu. S. Kivshar, Opt. Lett. 32, 3173 (2007). [27] D. Mihalache, D. Mazilu, F. Lederer and Yu. S. Kivshar, Opt. Lett. 32, 2091 (2007). [28] J. E. Heebner and R. W. Boyd, J. Mod. Opt. 49, 2629 (2002); P. Chak, J. E. Sipe and S. Pereira, Opt. Lett. 28, 1966 (2003). [29] J. J. Baumberg, P. G. Savvidis, R. M. Stevenson, A. I. Tartakovskii, M. S. Skolnick, D. M. Whittaker and J. S. Roberts, Phys. Rev. B 62, R16247 (2000); P. G. Savvidis and P. G. Lagoudakis, Semicond. Sci. Technol. 18, S311 (2003). [30] O. Morsch and M. Oberthaler, Rev. Mod. Phys. 78, 179 (2006). [31] I. Bloch, Nature Phys. 1, 23 (2005). [32] see the recent work of K. G. Makris and D. N. Christodoulides, Phys. Rev. E 73, 036616 (2006), as well as earlier ones such as K.Ø. Rasmussen, D. Cai, A.R. Bishop and N. Grønbech-Jensen, Phys. Rev. E 55, 6151 (1997) and P.G. Kevrekidis, I.G. Kevrekidis and B.A. Malomed, J. Phys. A 35, 267 (2002). [33] G. L. Alfimov, V. A. Brazhnyi, and V. V. Konotop, Physica D 194, 127 (2004). [34] M. I. Weinstein, Nonlinearity 12, 673 (1999). [35] D. E. Pelinovsky, P. G. Kevrekidis, and D. J. Frantzeskakis, Physica D 212, 20 (2005). [36] M. Golubitsky, D. G. Schaeffer, and I. Stewert, Singularities and Groups in Bifurcation Theory, Vol. 1 (Springer Verlag, New York, 1985). [37] D. E. Pelinovsky, P. G. Kevrekidis, and D. J. Frantzeskakis, Physica D 212, 1 (2005). [38] P. G. Kevrekidis, B. A. Malomed, D. J. Frantzeskakis and R. Carretero-Gonz´ alez, Phys. Rev. Lett. 93, 080403 (2004). [39] R. Carretero-Gonz´ alez, P. G. Kevrekidis, B. A. Malomed, and D. J. Frantzeskakis, Phys. Rev. Lett. 94, 203901 (2005). [40] M. Lukas, D. Pelinovsky and P. G. Kevrekidis, Physica D 237, 339 (2008).
