1. Introduction

Interplay between nonlinear effects and diffraction control in planar arrays of optical waveguides and in the light induced lattices in photorefractive materials have recently attracted substantial interest, see, e.g., [1] for a review. Spatiotemporal nonlinear effects relying on group velocity dispersion and, in particular, the topic of light bullets, is another active research area [2]. Despite numerous interesting effects, which interplay of controllable diffraction and GVD can bring into nonlinear optics, the experimental efforts devoted to these problems in recent years have been limited [3]. This is most likely because limitations imposed by the material and structure parameters. Photonic crystal fibers could be considered as a notable exception [4], which allows large degree of the GVD engineering by periodic variations of the refractive index, but on the other hand it most often enforces geometrical suppression of the diffraction, rather than allowing exploration of the functionality of the diffraction control as it has been done in planar structures [1]. On the theoretical side the results on the spatio-temporal and soliton effects in nonlinear waveguide arrays and related band-gap systems are not that numerous either. They include the existence of the continuous in time and discrete in space bullets [5], multi-dimensional solitons in systems with a gap-like dispersion along one coordinate and the free space diffraction along the other dimensions [6, 7], spatio-temporal X-waves in waveguide arrays [8] and solitons in coupled waveguides with Bragg gratings [9].

One of the theoretical methods of dealing with spatio-temporal problems, which has been well developed for the continuous systems, is to take a spatial soliton and consider its robustness with respect to the time dependent perturbations [10, 11, 12]. It is known that presence of the GVD usually destabilizes spatial solitons in the continuous systems via growth of the frequency side-bands. This process is usually called modulational instability (MI). Spatio-temporal dynamics induced by solitonic MI in continuous systems has not only been studied theoretically [10, 11, 12], but more recently observed experimentally [13, 14, 15]. Similar instabilities can be expected to occur for the discrete solitons in waveguide arrays. Our primary goal here is to consider temporal MI of the single peak spatial solitons in waveguides arrays with GVD. This problem becomes particularly relevant now because of the advances in fabrication of low loss planar silicon-on-insulator (SOI) structures with strong and controllable GVD and short coupling length, see, e.g., [16, 17, 18]. Most importantly the strong spatial confinement on the nano-scale and large Kerr nonlinearity bring the power levels required for the MI and soliton related effects to show on the millimeter to centimeter lengths down to few watts. The estimates for parameters typical for SOI waveguides show that all the effects described below are within the experimental reach, while GVD induced modulational instability in a single SOI waveguide has already been reported in [18].

We should also mention here that the MIs of spatially extended, i.e. non-localized, super-modes of the waveguide arrays induced by the discrete diffraction have been recently observed experimentally [19] and previously studied theoretically [20]. Papers [21, 22] studied the same case of diffraction induced MI, but for the supermodes consisting from the temporal solitons. Spatial MI of spatial surface solitons in optical lattices has been recently reported [23]. However, none of the known to us studies explored the problems of GVD induced instabilities of the spatial discrete solitons.

2. Model

We model an array of dielectric waveguides by a set of coupled NLS equations

with periodic boundary conditions U _N+1=U ₁ and U ₀=U_N . Here n enumerates the waveguides, τ and ζ are the time and coordinate along the waveguide, respectively. γ=2πn ₂/(Sλ) is the nonlinearity parameter, where S is the effective mode area and n ₂ is the Kerr coefficient. β ₂ is the GVD coefficient. κ=π/(2l_c ) is the coupling parameter and l_c is the coupling length. In order to put Eqs. (1) into dimensionless form we divide them by some fixed length l. Then in terms of dimensionless propagation distance z=ζ/l and dimensionless time $t=\frac{\tau }{\sqrt{\mid {\beta }_2\mid l}}$ Eqs. (1) take the form

$i{\partial }_z{A}_n-\frac{1}{2}s{\partial }_t^2{A}_n+C\left(A_{n+1}+A_{n-1}-2A_n\right)+{\mid A_n\mid }^2{A}_n=0,$ ,

with C=πl/(2l_c ), s=sign(β ₂), and $A_n=U_n\sqrt{\gamma l}.$ .

3. Stability analysis

Let us consider time-independent discrete soliton solutions of Eqs. (1) [20]. These solutions having the form U_n =a_ne^iqz can be found numerically, see Figs. 1(a,d,g). It is known that the discrete solitons are dynamically stable for β ₂=0 and the Vakhitov-Kolokolov stability criterion ∂_qQ>0, where $Q={\displaystyle {\sum }_{n=1}^N|a_n{|}^2}$, is satisfied for them [24]. As it is shown below the instabilities can arise as soon as GVD is taken into account.

To study stability with respect to time-dependent perturbations, we make the following ansatz

Here ω is the perturbation frequency. The linearized equations for the amplitudes of small perturbations ε _n,± can be transformed into the operator form

 

Fig. 1. MI of discrete solitons. The first row corresponds to C=7, the second to C=15 and the third one to C=30, respectively. q=10 and N=51 for all the panels. The right column shows transverse profiles of the discrete solitons. The middle column presents the frequency dependence of the MI growth rate (Iml>0) in the anomalous GVD regime (s<0). The right column shows all the unstable eigenvalues in the case of the normal GVD (s>0). Letters ‘N’ and ‘S’ mark the ‘neck’ and ‘snake’ instabilities, respectively.

Download Full Size | PDF

where ε⃗=(ε _1,+,ε _1,- …, ε _N,+,ε _N,-) ^T and L̂ ₁ is the diagonal N×N matrix: L̂ ₁=diag(1, -1, …, 1, -1). The matrix L̂ ₀ has the form

where q̃=q+2C. MI of discrete solitons manifests itself through a growth of perturbations in a certain range of frequencies ω. This means that in this range there exists an eigenvalue of the problem (3) such that Im(λ)>0. According to the classical results on MI of the bright solitons in the continuous NLS model [10] the eigenvectors of the corresponding eigenvalue problem can be either symmetric or antisymmetric with respect to the reflection about the soliton center. The instabilities associated with symmetric and antisymmetric eigenvectors are usually referred to as a ‘neck’-instability and ‘snake’-instability, respectively [10, 11, 12]. We will adopt the same terminology in our stability analysis of the discrete solitons. If the soliton is centered at n=n ₀, then the eigenvectors of Eq. (3) can be either symmetric on the replacement of n ₀+m with n ₀-m (as the soliton itself) or antisymmetric.

It follows from Eq. (3) that the stability of a discrete soliton with respect to the perturbations with ω=0 is determined by the eigenvalues of the matrix L̂ ₀. Since the soliton is stable in the absence of GVD, all the eigenvalues of L̂ ₀ are located on the real axis of the complex λ-plane. The zero eigenvalues of L̂ ₀, which are related to the continuous symmetries of Eqs. (1) can be used to get some analytical results. The symmetry transformation U_n →U_ne^iθ (where θ is an arbitrary phase) together with the Hamiltonian structure of these equations imply the existence of the two zero eigenvalues of the matrix L̂ ₀ with the same eigenvector x⃗ ₀. This means that the following two identities: L̂ ₀ x⃗ ₀=0 and L̂ ₀ x⃗ ₁=x⃗ ₀, which can be rewritten in the form ${\widehat{L}}_0^2$ x⃗ ₁=0, are satisfied. Here x⃗ ₀=(a ₁, -a ₁, …, a_N , -a_N ) ^T , x⃗ ₁=-∂_q (a ₁, a ₁, …, a_N , a_N ) ^T . When the frequency ω deviates from zero, see Eq. (3), these two eigenvalues move away from the origin in the complex λ-plane. They move in the opposite directions either along the real or along the imaginary axis. The latter case corresponds to the ‘neck’-instability associated with the vector x⃗ ₀ that has the same symmetry as the soliton itself. Note, that in the discrete case the continuous symmetry with respect to the lateral shifts of the soliton position is absent, and therefore the degeneracy of the zero eigenvalue of L̂ ₀ in the discrete model is half less than in the continuous one [10, 11, 12].

 

Fig. 2. The left column shows patterns of the ‘neck’ instability for anomalous GVD (s< 0) for 3 consequential values of the propagation distance z: C=7. The middle column shows patterns of the ‘neck’ instability for normal GVD (s>0): C=7. The right column shows patterns of the ‘snake’ instability for normal GVD (s>0): C=30. q=10 and N=51 for all the panels.

Download Full Size | PDF

Assuming that ω≪1 we write the following asymptotic expansions for the eigenvalue and corresponding eigenvector: λ=ωλ ₁+ω ² λ ₂ + … and x⃗=x⃗ ₀ +iλ ₁ ωx ₁ +ω ² x⃗ ₂ + …. Substituting these expansions into Eqs. (3), we get in the third order of the perturbation theory the solvability condition ${\lambda }_1^2$=sQ/(2∂_qQ). According to this condition the long wavelength instability of the ‘neck’-type takes place for s<0 (anomalous GVD). This result is not surprising because essentially the same instability persists even for κ=C=0. The middle column in Fig. 1 shows the ‘neck’ instability growth rates for different values of C. The low frequency part of the instability growth rates is approximately described by the above analytical expression. Direct numerical modeling of Eqs. (1) with the initial condition corresponding to a discrete soliton shows in this case formation of the regular trains of spatio-temporal quasi-solitons, which are discrete in space and continuous in time, see the left column in Fig. 2. Thus development of MI for the anomalous GVD case is qualitatively similar to the MI of the spatial solitons in the continuous NLS equation with saturable nonlinearity [25].

In the case of normal GVD (s>0) we have found complex instability spectra consisting from multiple sidebands, see the rightmost column in Fig. 1. These type of spectra appear to be specific to the discrete solitons. For the relatively small coupling strength, i.e. sufficiently far from the continuous limit, the dominant instability is of the ‘neck’-type. This instability leads to the break-up of the initial soliton to the localized lumps of light, which disperse with further propagation, see the middle column in Fig. 2. Contrary, for the continuous 2D NLS with the normal GVD the anti-symmetric ‘snake’-like instability dominates dynamics of 1D bright solitons [10, 26]. More, recent studies [26] have demonstrated that the symmetric ‘neck’-type eigenvectors also can be unstable in the continuous hyperbolic 2D NLS, with their growth rate been below the one for the ‘snake’ MI. As we have found in our model the ‘neck’ instability dominates the dynamics of the discrete solitons for small to moderate values of the coupling coefficient, see Fig. 1. Only for strong coupling, when the system becomes quasi-continuous, the ‘snake’ instability starts to be dominant over the ‘neck’ one. This leads to the break up of the discrete soliton in the snake-like fashion, see the right most column in Fig. 2. Note, that the dominant MI band of the ‘neck’-type found for small frequencies ω close to zero is associated with a pair of complex eigenvalues λ. At the same time the dominant band of the ‘snake’ MI and the strongest peak of the ‘neck’- MI at the relatively large frequencies ω have purely imaginary eigenvalues λ.

4. Physical estimates

As a guideline for physical estimates we consider parameters typical for SOI waveguides. In particular for a channel waveguide with width 480nm and thickness 220nm [16] GVD at 1.5µm is anomalous and its value is ≃ 580ps/nm/km. For width below 400nm or above 640nm GVD becomes normal. The coupling length for spacing around 400nm can be estimated at 200µm [17]. The kerr coefficient n ₂ for silicon is ~ 6×10^-14cm²/W. For our power estimates we take the effective area S≃0.3µm². Then one can show that C=7, 15, 30 used in numerical modeling give the following values for the dimensionless unit of z: 0.88mm, 2mm, 3.8mm. One unit of the time t corresponds to 22fs, 34fs, 47fs and unit of the peak power – to 1.7W, 0.8W, 0.4W, respectively. Considering that in the best SOI waveguides the loss is few dB/cm and remembering about two-photon and free carrier absorption, our power estimates should be scaled up. In particular, MI in a single SOI waveguide reported in [18] has been observed for 10W of pump power. Note, that though more detailed account of the above absorption mechanisms is desirable in future research, it can be forecasted, that they will simply proportionally suppress the instability growth rate, without qualitative changes in the effect itself. The above time estimates show that MI can be observed already with pico-second pump pulses, when the role of the free carrier absorbtion and dispersion is negligible.

5. Summary

We have analyzed modulational instability of bright discrete solitons in the waveguide arrays with group velocity dispersion. In the case of normal GVD we have found multiple instability bands. For weak to moderate strength of coupling the discrete solitons exhibit the ‘neck’ instability leading to breakup of the solitons into a train of dispersive pulses. Only for strong coupling, i.e., in the quasi-continuous limit, this instability is getting gradually suppressed by the ‘snake’ instability known for the 2D continuous NLS model [10]. In the case of anomalous GVD the expected neck type instability leads to formation of composite discrete-continuous spatio-temporal quasi-solitons.