Surface Solitons in Three Dimensions

Surface waves have been a subject of interest in a variety of contexts, including surface plasmons in conductors Barnes and optical solitons in waveguide arrays Christo in physics, surface waves in isotropic magnetic gels Bohlius 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 Makris and observed experimentally Suntsov that discrete localized nonlinear waves can be supported at the edge of a semi-infinite array of nonlinear optical waveguide arrays. Such solitary waves were predicted to exist not only in self-focusing media (as in the above-mentioned works), but also between uniform and self-defocusing media Makris ; kartprl , or between self-focusing and self-defocusing media (e.g. in debra ). They have been subsequently observed in media with quadratic Sivilogou and photorefractive rosberg ; smirnov nonlinearities. In the two-dimensional (2D) geometry, stable topological solitons have been predicted in a saturable medium Kartashov , which constitute generalizations to lattice vortex solitons predicted in Ref. Malomed . Quasi-discrete vortex solitons have been experimentally observed in a self-focusing bulk photorefractive medium Neshev . Theoretical predictions for a variety of species of discrete 2D surface solitons Susanto ; Makris2 ; Kartashov2 ; BluKon ; Vicencio , corner modes Makris2 ; BluKon , as well as surface breathers BluKon , were reported. Subsequent work has resulted in the experimental observation of 2D surface solitons, both fundamental ones and multi-pulse states, in photorefractive media christo2d , as well as in asymmetric waveguide arrays written in fused silica kart2d . Recently, surface solitons in more complex settings, such as chirped optical lattices in 1D and 2D kart_mol1 ; kart_mol2 , at interfaces between photonic crystals and metamaterials shadrivov , and in the case of nonlocal nonlinearity moti1 ; moti2 , 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. dumitru (the respective 2D surface structures were examined in Ref. dumitru1 ).

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 earlier1 , polaritons earlier2 , or Bose-Einstein condensates in the vicinity of an edge of a strong 3D optical lattice morsch ; bloch . 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. Susanto . 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 structure. 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. Christodoulides .

To formulate the model, we introduce unit vectors $\mathbf{e}_{1}=(1,0,0)$, $\mathbf{e}_{2}=(0,1,0)$, and $\mathbf{e}_{3}=(0,0,1)$ and define lattice sites by $\mathbf{n}=\sum_{j=1}^{3}n_{j}\mathbf{e}_{j}$ with integer $n_{j}$. We assume that the lattice occupies a semi-infinite space, $n_{3}\geq 1$, and its dynamics obeys the discrete nonlinear Schrödinger (DNLS) equation in its usual form,

Here $\phi_{\mathbf{n}}$ is a complex discrete field, $\varepsilon$ is the coupling constant, $\dot{\phi}_{\mathbf{n}}$ stands for the time derivative, the parameter $\sigma=\pm 1$ determines the sign of the nonlinearity (focusing or defocusing respectively), and $\Delta\phi_{\mathbf{n}}$ is the 3D discrete Laplacian:

for $n_{3}\geq 2$, while for $n_{3}=1$ the term with subscript index $\mathbf{n}-\mathbf{e}_{3}$ is to be dropped (note that $\mathbf{e}_{3}$ 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 $n_{3}$ up to $n_{3}=-\infty$, supplementing the equation with the anti-symmetry condition,

Indeed, this condition implies $\phi_{n_{1},n_{2},0}\equiv 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, $\phi_{\mathbf{n}}\rightarrow 0$ at $n_{1,2}\rightarrow\pm\infty$ and $n_{3}\rightarrow\infty$. Additionally, $s=\pm 1$ in Eq. (2) —this parameter is introduced for convenience (see Sec. III.2) and can be freely rescaled using the transformation $\phi\rightarrow\phi\,e^{i\nu t}$ for an appropriate choice of $\nu$ and time rescaling. Stationary solutions to Eq. (1) will be sought for as $\phi_{\mathbf{n}}=\exp{(i\Lambda t)}u_{\mathbf{n}}$, where $\Lambda$ is the frequency and the lattice field $u_{\mathbf{n}}$ obeys the equation

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 $\varepsilon$, 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.

First, we outline some general properties of the model. Equation (1) conserves two dynamical invariants, namely the norm $N$,

and the Hamiltonian $H$,

where the asterisk stands for complex conjugation. Stationary solutions to Eq. (4) with $\sigma=\pm 1$ are connected by the staggering transformation ABK ; BluKon : if $u_{\mathbf{n}}$ is a solution for some $\Lambda$ and $\sigma=+1$, then $(-1)^{n_{1}+n_{2}+n_{3}}u_{\mathbf{n}}$ is a solution for $\widetilde{\Lambda}=12s-\Lambda$ and $\sigma=-1$. Consequently, it is sufficient to perform the analysis of stationary solutions, including their stability, for a single sign of the nonlinearity; thus, below we will fix $\sigma=+1$ (corresponding to the case of onsite self-attraction).

Solutions to Eq. (4) in half-space $n_{3}\geq 1$, subject to boundary condition $\phi_{\mathbf{n}}=0$ for $n_{3}=0$, as defined above, may be continued anti-symmetrically for the entire 3D space by setting $U_{\mathbf{n}}\equiv u_{\mathbf{n}}$ for $n_{3}\geq 1$ and $U_{\mathbf{n}}\equiv-u_{\mathbf{n}}$ for $n_{3}\leq-1$. Then, according to results of Ref. Weinstein , this leads to an immediate conclusion, namely that there exists a minimum norm $N_{\min}$ necessary for the existence of localized surface states in the present model. In other words, no surface modes survive in the limit of $N\rightarrow 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 BluKon . 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 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 $\varepsilon$ 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.2).

To find solution families, we start with the anti-continuum (AC) limit, $\varepsilon=0$ Pelinovsky . 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 $\varepsilon>0$ is determined by the Lyapunov’s reduction theorem Golubitsky . More specifically, the solution is expanded as a power series in $\varepsilon$, 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 Pelinovsky .

The linear stability is then studied, starting from the usual form of the perturbed solution,

where $\delta$ is a formal small parameter, and $\lambda$ is a stability eigenvalue associated with eigenvector $\psi=\{a_{\mathbf{n}},b_{\mathbf{n}}^{\ast}\}$. Substituting this into Eq. (1) yields the linearized system,

which can be cast in the form

where $\mathcal{A}$ and $\mathcal{B}$ are vectors composed by elements $a_{\mathbf{n}}$ and $b_{\mathbf{n}}^{\ast}$, respectively, while the matrices $\mathcal{H}^{(p,q)}$ ($p,q\in\{1,2\}$) are given by,

An underlying stationary solution is (spectrally) unstable if there exists a solution to Eq. (8) with $\mathrm{Re}(\lambda)>0$. Otherwise, the stationary solution is classified as a spectrally stable one. As explained in Ref. Pelinovsky1d , the Jacobian of the above mentioned solvability conditions is intimately connected to the full eigenvalue problem. More specifically, if the eigenvalues $\gamma$ of the $M\times 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 $\lambda=\sqrt{2\gamma}\varepsilon^{p/2}$, where $p$ is the number of lattice sites that separate the adjacent excited nodes of the configuration at the AC limit.

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 $\varepsilon$ 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 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 $\lambda$ 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 $\Lambda$, only at values of the coupling constant $\varepsilon\approx\Lambda$ earlier3 . Another structure with a relatively wide stability region is the dipole (the more stable the wider the distance between its constituent sites vor3dnew ), 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 values of $\varepsilon$, 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 $\varepsilon>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 $\varepsilon<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 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 $\sim 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.

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 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. Pelinovsky1d ; Pelinovsky ; Pelinovsky3d for the analysis of localized states in 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 higher-order 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.