Embedded solitons in second-harmonic-generating lattices
aa r X i v : . [ n li n . PS ] F e b Embedded solitons in second-harmonic-generating lattices
H. Susanto
Department of Mathematics, Khalifa University, Abu Dhabi Campus, PO Box 127788,Abu Dhabi, United Arab EmiratesDepartment of Mathematical Sciences, University of Essex, Wivenhoe Park,Colchester CO4 3SQ, United Kingdom
Boris A. Malomed
Department of Physical Electronics, School of Electrical Engineering, Faculty of Engineering, andCenter for Light-Matter Interaction, Tel Aviv University, Tel Aviv 69978, IsraelInstituto de Alta Investigaci´on, Universidad de Tarapac´a, Casilla 7D, Arica, Chile
Abstract
Embedded solitons are exceptional modes in nonlinear-wave systems with the propaga-tion constant falling in the system’s propagation band. An especially challenging topicis seeking for such modes in nonlinear dynamical lattices (discrete systems). We addressthis problem for a system of coupled discrete equations modeling the light propagation inan array of tunnel-coupled waveguides with a combination of intrinsic quadratic (second-harmonic-generating) and cubic nonlinearities. Solutions for discrete embedded solitons(DESs) are constructed by means of two analytical approximations, adjusted, severally,for broad and narrow DESs, and in a systematic form with the help of numerical calcu-lations. DESs of several types, including ones with identical and opposite signs of theirfundamental-frequency and second-harmonic components, are produced. In the mostrelevant case of narrow DESs, the analytical approximation produces very accurate re-sults, in comparison with the numerical findings. In this case, the DES branch extendsfrom the propagation band into a semi-infinite gap as a family of regular discrete soli-tons. The study of stability of DESs demonstrates that, in addition to ones featuring thewell-known property of semi-stability, there are linearly stable DESs which are genuinelyrobust modes.
Keywords: embedded soliton, discrete soliton, variational approximation
1. Introduction
One of basic scenarios of nonlinear optics is the generation of second-harmonic (SH)waves by the fundamental-frequency (FF) input in waveguides with a quadratic, alias χ (2) , material response [1]. The interplay of the χ (2) nonlinearity with linear paraxial Email address: [email protected] (H. Susanto)
Preprint submitted to Chaos, Solitons, & Fractals February 17, 2021 iffraction of light in planar waveguides gives rise to spatial solitons of the χ (2) type,which are built of quadratically coupled FF and SH components. Such solitons werefirst predicted as particular exact solutions of a system of two coupled equations forthe paraxial propagation of amplitudes of the FF and SH waves [2], and later createdexperimentally [3]. Theoretical and experimental studies of χ (2) solitons have grown intoa vast research area [4, 5, 6].Varieties of χ (2) dynamical effects and, in particular, soliton modes, are additionallyexpanded in settings which combine quadratic and cubic ( χ (3) ) terms [7, 8]. There are twodistinct types of such combinations: competing, if the χ (3) term is self-defocusing, andcooperating, if the latter term has the focusing sign. One of specific phenomena predictedin models with the combined χ (2) : χ (3) nonlinearity is the existence of embedded solitons [9] i.e., localized modes whose propagation constant, rather than belonging, as usual, tospectral gaps, in which the propagation of linear waves is suppressed by the system’sdispersion relation, falls (is embedded ) in a propagation band. Normally, solitons cannotexist in such a case, as the resonance between propagation constants of the solitons andlinear waves gives rise to decay of the (quasi-) soliton modes into radiation. Nevertheless,in exceptional cases. i.e., at some isolated values of the soliton’s propagation constant,the decay rate may vanish (in other words, the resonance strength becomes exactlyequal to zero), which explains the existence of embedded solitons. This possibility wasoriginally proposed in Refs. [10] and [11] (without naming such modes “embedded” ones),see an early review of the topic in Ref. [12]. In the linear version of the χ (2) system,dispersion relations for the FF and SH waves are decoupled, each consisting of a semi-infinite propagation band and forbidden gap. Boundaries between the gap and band inthe FF and SH spectra are mutually shifted, in the general case, due to the mismatchbetween the FF and SH components. In terms of the FF wave, the propagation constantof any soliton must always belong to the respective gap. However, the SH propagationequation, with the SH amplitude driven by the square of the FF, may have nonlinearizable solutions. For such states, the quadratic (FF × FF) term in the SH equation cannot beomitted when considering exponentially decaying tails of the soliton’s fields, as it hasthe same order of magnitude as linear SH terms (“tail-locked” modes) [9]. Thus, the SHpropagation constant of these solutions does not need to obey the spectral structure ofthe linearized equation, falling, counter-intuitively, into the propagation band for the SHcomponent.Because, as mentioned above, embedded solitons exists only at specially selectedvalue(s) of the propagation constant, they do not form continuous families, being cate-gorized, accordingly, as codimension-one solutions (see, e.g. Ref. [13]). Another specificfeature of embedded solitons is their (in)stability. In most cases, they are neutrally sta-ble against small perturbations in the linear approximation, but are unstable against the(slower-than-exponential) growth of perturbations driven by terms which are quadratic,with respect to the small perturbations [9]. In some cases, such “semi-stable” solitonsmay seem as practically stable ones, if the nonlinear growth of the perturbations isextremely slow [12], and there are particular models which support completely stableembedded solitons [14].One of ramifications of studies of χ (2) solitons is the creation of their discrete coun-terparts in arrays (lattices) of tunnel-coupled waveguides with the quadratic or combinedquadratic-cubic intrinsic nonlinearities in individual guiding cores [15, 16]. A challengingpossibility is to search for discrete embedded solitons (DESs) in the respective system2f lattice-propagation equations for amplitudes of the FF and SH waves, coupled by thequadratic nonlinearity acting at lattice sites. This objective was a subject of Ref. [17].In that work, only a partial consideration of the problem was performed; in particular,only DESs with opposite signs of the FF and SH components were found. Note that noanalytical methods were applied in Ref. [17]. In the present work, we aim to developthe consideration of DESs in χ (2) lattices, elaborating an analytical approximation (theaveraging method), in parallel to producing systematic numerical results, and producingDESs of more general types, featuring both identical and opposite signs of the FF andSH components. Additionally, we also consider stability of the DES.The model is formulated in Section 2, and two analytical methods are presentedin Section 3. One of them (the quasi-continuum approximation) is relevant for broadDESs, and the opposite anti-continuum approximation, applies, with high accuracy, tostrongly discrete solitons. Numerical results are reported, and compared to the analyticalapproximations, in Section 4. Stability of DES is investigated in the same section bymeans of computing eigenvalues of the linearisation operator and direct simulations ofthe perturbed evolution. We also report the existence of linearly stable embedded solitonsin the weakly coupled lattice. The paper is completed by Section 5.
2. The model
We start with the discrete one-dimensional system with the on-site χ (2) : χ (3) nonlin-earity, which was introduced in Ref. [17]: i du n dz + ǫ ∆ u n + u ∗ n v n + γ ( | u n | + 2 | v n | ) u n = 0 , (1a) i dv n dz − ǫδ ∆ v n + qv n + 12 u n + 2 γ ( | v n | + 2 | u n | ) v n = 0 . (1b)Here ǫ is the coefficient of coupling between adjacent sites of the underlying lattice forcomplex amplitude u n of the FF wave, δ ≷ v n of the SH wave, the finite-difference Laplacian is ∆ X n = X n +1 + X n − − X n , the SH-FF mismatch is accounted for by constant q , the coefficient of the χ (2) interaction is scaled to be 1, and γ , are coefficients of the χ (3) self-interaction for theFF and SH waves. All coefficients in Eqs. (1a) and (1b) are real, while ∗ stands for thecomplex conjugate.Note that, rescaling the lattice fields u n and v n by a common factor, one can fix thevalue of γ , while γ remains an independent parameter. In this work, we use this optionby fixing | γ | = 0 .
05. As for γ , we produce numerical results for γ = γ , as variationof ratio γ /γ does not lead to a conspicuous change of the results. The sign of γ , is essential, positive and negative ones representing, respectively, the self-focusing anddefocusing χ (3) nonlinearity. In the absence of the χ (2) terms, this sign may be fixedby means of the staggering transformation , ( u n , v n ) → ( − n ( u n , v n ) [20]. However,the present system does not admit it, as the sign of the χ (2) coefficient is already fixed.Below, we consider both signs of γ , .The FF lattice coupling constant in Eq. (1a) may be set to be positive, ǫ > δ > δ < u n = e iωz U n , v n = e iωz V n , (2)where real ω is a propagation constant, and real discrete fields U n , V n satisfy the followingequations: ( − ω + ǫ ∆) U n + U n V n + γ (cid:0) U n + 2 V n (cid:1) U n = 0 , (3a)( q − ω − ǫδ ∆) V n + 12 U n + 2 γ (cid:0) V n + 2 U n (cid:1) V n = 0 . (3b)Llinearizing these equations in a formally complex form, it is straightforward to derivedispersion relations for the plane waves ∼ exp ( ikn ) in the FF and SH components: ω = − ǫ (1 − cos k ) , ω = ǫδ (1 − cos k ) + q/ , (4)where k is the spatial wavenumber. As it follows from Eq. (4), the propagation band forthe plane waves in the FF component is − ǫ < ω < , (5)while for the SH component it is q/ < ω < q/ ǫδ for δ > , (6) q/ − | δ | ǫ < ω < q/ δ < . (7)Spectral intervals not covered by the propagation bands are considered as gaps. Thus,each band given by Eqs. (5) and (6), (7) is located between two semi-infinite gaps forthe FF and SH components, respectively. The gaps extending to ω → + ∞ and −∞ mayhost, severally, discrete solitons of the unstaggered and staggered types, the latter oneassuming the structure of the lattice fields ∼ ( − n [20]. To look for DES solutions, weassume that the FF propagation constant belongs to the semi-infinite gap correspondingto unstaggered discrete solitons, i.e., ω >
0, as it follows from Eq. (5). Simultaneously, itis assumed that the SH propagation constant, 2 ω , falls into the corresponding propagationband given by Eqs. (6) or (7) (otherwise, one would be dealing with regular discretesolitons, rather than DESs).
3. Analytical approximations ǫ ≫ ) For ǫ ≫
1, equations (3a) and (3b) may be approximated by their continuous coun-terparts: − ωU + U xx + U V + γ (cid:0) U + 2 V (cid:1) U = − ǫ U xxxx + O (1 /ǫ ) , (8a) − ωV − δV xx + qV + 12 U + 2 γ (cid:0) V + 2 U (cid:1) V n = − ǫ V xxxx + O (1 /ǫ ) . (8b)4n this limit (if the fourth derivatives may be omitted), and under the assumption thatthe SH component is much weaker than the FF one, i.e., | V | ≪ | U | Eqs. (3a), (3b)give rise to a known single-humped embedded-soliton solution [9].The variational approximation, which may be efficiently used to construct soliton-like solutions of discrete equations [22, 23, 24], does not apply to Eqs. (3a), (3b), aswell as to their continuous-limit counterparts (8a), (8b), because these systems cannotbe derived from Lagrangians, except for the case of γ = 2 γ (numerical results areproduced below for the case of γ = γ , when the Lagrangian representation is notavailable). Nevertheless, it is possible to apply the averaging method [21], which uses anadopted ansatz, such as the one introduced below in Eq. (9), as an “optimal” approximatelocalized solution of Eqs. (3a), (3b). The averaging method amounts to the variationalone if the underlying equation is Lagrangian [21].Motivated by the example of the exact embedded soliton found in Ref. [9], we assumethat localized solutions of Eqs. (8a) and (8b) may be approximated by U = A sech (cid:0) √ ωx (cid:1) , V = B sech (cid:0) √ ωx (cid:1) . (9)Following the averaging method, we substitute the ansatz in Eqs. (8a), (8b), multiplythe first equation by ∂U/∂A and the second one by ∂V /∂B , and integrate the resultingexpressions over −∞ < x < + ∞ , to derive the following relations between parametersof the ansatz: A = 1 γ (cid:18) ω − B γ − B − ω ǫ (cid:19) , (10a) A = 2 B (cid:2)
42 (2 δ − ω − B γ − q − ǫ − ω (cid:3) Bγ + 5) . (10b)Approximate solutions given by Eqs. (10a), (10b) are candidates for embedded solitonsin the case of δ >
0. However, they do not select true DESs, as the actual SH componentwould likely include tails which do not vanish at x → ±∞ . We therefore need to imposean additional condition, which is the orthogonality of the localized solution (9) to thetail [25, 26].Because of the possible presence of the nonvanishing tail in the SH component, weseek for solutions in the form of V = V + CV , (11)where V is approximated by Eq. (9), C is a constant, and the asymptotic form of thetail has the form of V ∼ e ± ikx at x → ±∞ , where wavenumber k is given by dispersionrelation following from the linearization of Eq. (8b), i.e., k = 6 ǫ (cid:16) − δ ± p δ + (2 ω − q ) / ǫ (cid:17) . (12)Making use of the averaging method once again, we multiply Eq. (8b) by ∂V /∂C , sub-stitute C = 0 in the resulting expression (as we want to identify a DES), and integratethe equation over −∞ < x < + ∞ , to derive the orthogonality condition J ≡ Z + ∞−∞ N ( U , V ) cos( kx ) dx = 0 , (13)5here U is ansatz (9) for the approximate localized solution, and N ( U, V ) = − ωV − δV xx + qV + (1 / U + 2 γ (cid:0) V + 2 U (cid:1) V n + (12 ǫ ) − V xxxx . (14)The integral in Eq. (13) can be evaluated by means of residues of the analytically con-tinued integrand, Z + ∞−∞ N ( V ) cos( kx ) dx = − π X Im Res (cid:2) N ( V ) e ikx (cid:3) , (15)where the summation is to be performed over all poles that lie in the upper complexhalf-plane of x . The analytical function written in Eq. (9) has infinitely many polesat imaginary values of the integration variable, x = (1 + 2 n ) iπ/ √ ω , n = 0 , , , . . . .However, the asymptotic limit is determined solely by the contribution from the lowestpole with n = 0, which yields J ≈ − kπ ω exp (cid:18) kπ √ ω (cid:19) × (cid:20) γ
20 ( k + 4 ω )( k + 16 ω ) B + 32 A ω +2 ω B (cid:18) − ω + (cid:18) A γ + 32 δk + 32 q (cid:19) + 1 ω A k γ + 18 ǫ k (cid:19)(cid:21) . (16)Thus, Eqs. (10a), (10b), (12) and (16) predict values of parameters for DESs in thecontinuum limit, ǫ ≫ ǫ ≪ ) In the anti-continuum limit, weakly coupled solutions of Eqs. (3a), (3b) may beapproximated by the simplest ansatz [22], U n = A exp ( − α | n | ) , V n = B exp ( − α | n | ) , (17)where amplitudes A and B are free parameters, and α is determined by the conditionthat the ansatz must agree with a solution of the linearized version of Eq. (3a) for thetail of the FF component [25], i.e., with the first equation of system (4), in which α = ik is substituted: α = arccosh (1 + ω/ ǫ ) . (18)Next, we determine values of A and B in ansatz (17), again using the averaging method:the ansatz is substituted in Eqs. (3a), (3b), multiplying the first equation by ∂U n /∂A ,the second one by ∂V n /∂B , and performing the summation over −∞ < n < + ∞ , toobtain (cid:0) A γ + B (cid:1) F + B γ F + ǫe − α − ω/ − ǫ = 0 , (19a) (cid:18) q − ω + A B (cid:19) F + 2 B γ F + 4 A γ F + 2 ǫδF = 0 , (19b)where F n ≡ coth ( nα/ B may also be approximated by demanding that, if theansatz is inserted in Eq. (3b), the tail of the SH component must be “locked” to the6quared FF tail at | n | → ∞ (i.e., the tails of both components must jointly solve theasymptotic form of the equation): B ≈ A [ ω − q/ ǫδ (cosh(2 α ) − − . (20)However, the latter approximation overestimates the role of the DES’ tails in comparisonto its core region, which limits the applicability of the approximation, therefore we donot employ it here.Similar to the situation considered above, Eqs. (19a) and (19b) do not select trueDESs, while the appropriate selection is provided by the condition of the orthogonalityof the soliton ansatz to the nonvanishing tail. To this end, the ansatz including the tail iswritten as V n = B exp ( − α | n | ) + C e V n , where C is a constant, e V n ∼ cos( nk ) as n → ±∞ ,and k is determined by the second equation of the system of dispersion relations 4, cf.Eq. (11). Also similar to what was done above, the orthogonality condition is obtainedby multiplying expression (3b) by e V n and performing the summation over n : J = + ∞ X n = −∞ cos ( kn ) × (cid:20) ( q − ω − ǫδ ∆) V n + 12 U n + 2 γ (cid:0) V n + 2 U n (cid:1) V n (cid:21) = 0 . (21)The infinite sum in Eq. (21) can be calculated by means of the following formula: + ∞ X n = −∞ cos ( kn ) e − jα | n | ≡ Re ( + ∞ X −∞ e − jα | n | + ik | n | ) = G j , (22)where G j = sinh ( jα ) / (cid:2) cosh ( jα ) − cos( k ) e − jα (cid:3) . Thus, Eq. (21) gives rise to the follow-ing selection condition for DES: J = (cid:2) Bδǫ (1 − cos k ) − ωB + qB + A / (cid:3) G + 2 B γ G + 4 A Bγ G = 0 . (23)The form of Eq. (23) clearly demonstrates that that this condition may be satisfied only in the case of the combined χ (2) : χ (3) nonlinearity: if the cubic terms are absent, i.e., γ = 0, Eq. (23) cannot hold.Thus, Eqs. (19a), (19b) and (23) predict parameter values for DESs in the anti-continuum limit, | ǫ | ≪
4. Numerical results
The numerical solution of stationary equations (3a) and (3b) was obtained by employ-ing an iterative optimization algorithm for the solution of nonlinear least-squares prob-lems through the function fsolve of Matlab . Once DES components U n and V n hadbeen found, we studies stability of the solution by substituting a perturbed solution, u n =( U n +(ˆ r n + i ˜ r n ) e λz +(ˆ r ∗ n + i ˜ r ∗ n ) e λ ∗ z ) e iωz and v n = ( V n +(ˆ s n + i ˜ s n ) e λz +(ˆ s ∗ n + ˜ s ∗ n ) e iλ ∗ z ) e iωz into Eqs. (1a), (1b) and linearizing them with respect to small perturbations ˆ r n , ˜ r n , ˆ s n and ˜ s n , arriving at the following problem for stability eigenvalue λ : − λ D ˆ r n ˜ r n ˆ s n ˜ s n = L , − U n L , + U n + 4 γ U n V n U n L , − U n + 8 γ U n V n L , + ˆ r n ˜ r n ˆ s n ˜ s n , (24a)7here we have defined a diagonal matrix,D ≡ − − , (25)and linearisation operators L , ∓ = ǫ ∆ − ω ∓ V n + γ (cid:0) (2 ∓ U n + 2 V n (cid:1) , (26) L , ∓ = − ǫδ ∆ − ω + q + 2 γ (cid:0) (2 ∓ V n + 2 U n (cid:1) . (27)DES is unstable when there are eigenvalues λ with nonzero real parts, and linearly stableotherwise. A semi-analytical approximation for identifying critical eigenvalues by meansof the variational formulation is possible (see, e.g., Ref. [27] and references therein), butwe do not pursue this direction here.After identifying the (in)stability of DES in terms of the eigenvalues, their perturbedevolution was investigated, simulating Eqs. (1a) and (1b) by means of the fourth-orderRunge-Kutta method in the temporal direction. First, in Fig. 1(a) we display an example of a numerically found DES, obtained in thelattice built of N = 100 sites, and its counterpart. produced by means of the approxima-tion based on Eqs. (9), (10a), (10b), (13), and (12), (16), for values of parameters in Eqs.(1a), (1b) ǫ = 10, q = 1, γ = γ = − .
05 (recall that γ , < χ (3) nonlinearity) and propagation constant ω = 0 .
7. Because DES is a solu-tion of the codimension-one type, for these arbitrarily chosen parameters the numericalsolution can be found only for a specially selected value of the remaining one, viz ., therelative strength of the lattice coupling for the SH wave, δ ≈ . δ ≈ . δ > . < ω < . ω = 0 . ǫ = 10 is large enough for this purpose. Note also that the relatively small valueof the propagation constant, ω = 0 .
7, makes the DES broad enough [as per Eq. (9)],which additionally justifies the use of the quasi-continuum formalism in this case. It isremarkable that the analytical approximation produces a virtually exact DES solutionin Fig. 1(a).Fixing parameters q = 1 and γ , = − .
05, and propagation constant ω = 0 . ǫ makes it possible to produce two branches of numerical solutionsof Eqs. (3a), (3b) in the ( ǫ, δ ) plane, as shown in Fig. 1(b) by the solid curve with aturning point at ǫ = ǫ min ≈ .
14, at which the two branches merge. Thus, Fig. 1(b)shows that the continuation of the DES families does not exist if the coupling constantfalls below the minimum value. Note that, for the currently fixed values, q = 1 and ω = 0 .
7, the SH band, given by Eq. (6), amounts to ǫδ > .
1. Obviously, the entire8 n U n V n (a) (b)
30 35 40 45 50 55 60 65 70 n V n U n (c)Figure 1: (a) A discrete embedded soliton (DES) for parameters of Eqs. ( 1a) and (1b) ǫ = 10, q = 1, ω = 0 . γ = γ = − . δ ≈ . ω = 0 . δ ≈ . ǫ, δ ), along which the numericalsolution produces DES states, with other parameters fixed at the same values as in panel (a). Thedashed curve shows a similar result produced by the analytical approximation. (c) A numerically foundDES solution belonging to the lower branch in panel (b), for the the same parameters as in (a), exceptfor δ ≈ . U n and V n , respectively. DES families represented in this figure are completely unstable. see Fig. 3 below. δ > δ min ≈ .
36 [in particular, DESs do not exist at δ <
0, when the couplingconstants in the FF and SH sublattices have the same signs, see Eqs. (1a) and (1b)].The analytical approximation based on Eqs. (9), (10a), (10b), (13), and (16) producesa single existence branch, shown by the dashed line in Fig. 1(b), which cannot beextended beyond the left end point, i.e., the approximation predicts, in a qualitativelycorrect form, only one branch of the DES solutions. We assume that this happens becauseansatz (9) is not sifficiently versatile to capture the other branch of the solutions. A betteransatz may be tried, but it would lead to an extremely cumbersome analysis.
10 12 14 16 18 20 n -5051015202530 U n V n (a) (b)Figure 2: The same as Fig. 1(a,b), but for the weakly coupled lattice with parameters q = 60, γ = γ =0 .
05, and propagation constant ω = 29. In panel (a), ǫ = 5 and δ ≈ − . U n and V n , are positive and negative ones. This DES is a completely stable one, see Fig. 5 below. Theshaded region in panel (b) shows the domain of existence of regular (non-embedded) solitons. The above results were obtained for a situation close to the continuum limit. Figure2 presents typical findings for the opposite case, close to the anti-continuum limit, in thelattice composed of N = 30 sites. First, an example of a strongly discrete embeddedsoliton is displayed in panel (a) for parameters of Eqs. (1a) and (1b) chosen as ǫ = 5, q = 60, γ = γ = 0 .
05, and propagation constant ω = 29, while the value of coefficient δ , at which the nongeneric DES is found numerically, is δ ≈ − . δ implies identical signs of the inter-site coupling in the FF and SHsublattices, on the contrary to the case displayed above in Fig. 1. At these values of theparameters, the SH band, defined by Eq. (7), is 24 . < ω <
30, to which ω = 29 indeedbelongs, sitting close to its edge.The analytical approximation for the strongly discrete DES, based on Eqs. (17), (18),and (23) is also displayed, for the same parameters, in Fig. 2(a). It is seen to be virtuallyidentical to the numerical counterpart.In a still simpler approximation, which may be applied to a strongly discrete soliton,its FF and SH amplitudes, U and V , can be found by completely neglecting the coupling10f the central site, n = 0, to the adjacent ones, n = ± U ≫ V , the following values: V ≈ − / (8 γ ) = − . , U ≈ p ( ω − V ) /γ ≈ . , (28)which are close enough to ones observed in Fig. 2(a).A family of DES modes, produced by the numerical solution of Eqs. (3a), (3b) for thesame parameters as in Fig. 2(a), i.e., q = 60, γ , = 0 .
05, and ω = 29, but with a varyingcoupling constant ǫ , and δ adjusted accordingly, is displayed in Fig. 2(b), along with itscounterpart produced by the above-mentioned analytical approximation, based on Eqs.(17), (18), and (23), It is clearly seen that the approximation is very accurate in thiscase, even for relatively large values of ǫ , for which the applicability of the anti-continuumformalism is not evident. Unlike the family shown above in Fig. 1, with identical signsof the FF and SH components, the present one always features opposite signs.Lastly, we note that, for the parameters chosen in Fig. 2(b), Eq. (7) defines the SHband as the area with | δ | ǫ > .
5. This condition means that the family presented in thisfigure belongs to the band at ǫ > ǫ ≈ . ǫ < ǫ immersesinto the bandgap, as a family of regular solitons (non-embedded ones). In the figure, theexistence region of the regular solitons is shaded. We address the stability of DESs by numerically solving the eigenvalue problem (24a).First, we consider the case of opposite signs of the coupling constant in the SH and FFlattices, i.e., δ > -0.05 0 0.05
Re( ) -50050 I m () (a) -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Re( ) -3 -50050 I m () (b)Figure 3: Panels (a) and (b) show spectra of complex stability eigenvalues for DES displayed in Figs. 1aand 1c, respectively. Having established the instability of the DESs shown in Fig. 1a, we display sim-ulations of their dynamics, under the action of two different perturbations, in Fig. 4.11irst, Fig. 4a demonstrates that the soliton, in spite of the presence of many unstableeigenvalues in its perturbation spectrum, as seen in Fig. 3a, is actually quite robust indirect simulations which include small random perturbations in the initial conditions. Inthis connection, we note that, with ǫ = 10 and δ close to 1 (the values used in Fig. 1),and the soliton’s width W ∼
10 (as seen in the same figure), the characteristic diffraction(Rayleigh) length is Z diffr ∼ W /ǫ ∼
10, i.e., the total propagation distance displayed inFig. 4(a) is ∼
30 diffraction lengths, which is actually a very large value in a real experi-ment [3]-[6]. On the other hand, if the initial state is perturbed merely by multiplying itby factor 0 .
99, the DES originally stays nearly intact, over the distance Z ≃ ∼ Z diffr ,and then quickly decays, as shown in Fig. 4b. The latter dynamical scenario is typical forthe “one-sided” semistability of embedded solitons [9], initiated by slow sub-exponentialgrowth of perturbations, which eventually leads to fast destructiom of the solitons. Theperturbed evolution of the DES displayed in Fig. 1c exhibits similar dynamical behavior(not shown here in detail). |u n |
30 40 50 60 70100200300 z |v n |
30 40 50 60 70 n (a) |u n |
30 40 50 60 70204060 z |v n |
30 40 50 60 70 n (b)Figure 4: The evolution of DES from Fig. 1a(a) under the action of complex random-noise perturbationswith maximum modulus amplitude 2 × −
2; and (b) after being initially perturbed by rescaling itsshape, according to u n ( t = 0) = 0 . × U n , v n ( t = 0) = 0 . × V n . Next, we consider the case of δ <
0, i.e., opposite signs of the intersite linear couplingsin the FF and SH lattices in Eqs. (1a) and (1b). For the DES shown in Fig. 2a, its linear-perturbation spectrum is plotted in Fig. 5a, which does not include unstable eigenvalues(the spectrum is purely imaginary). This is the case with all DESs along the existencecurve plotted in Fig. 2b. Direct simulations of the evolution of the soliton subjected tothe scaling perturbation, similar to that which led to the destruction of the semistablesoliton in Fig. 4(b), but actually much stronger (with the initial scaling factor 0 . . δ < Re( ) -60-40-200204060 I m () (a) |u n | z |v n | n (b)Figure 5: (a) The spectrum of perturbation eigenvalues for the DES from Fig. 2a. (b) The simulatedevolution of the same soliton, with an initial perturbation in the form of u n ( t = 0) = 0 . × U n , v n ( t =0) = 0 . × V n . This, relatively strong, perturbation leaves the soliton completely stable, unlike a muchweaker initial scaling perturbation which destroys the DES in Fig. 4(b).
5. Conclusion
In this work, we have revisited the problem of the creation of DESs (discrete em-bedded solitons) in the model of the light propagation in an array of tunnel-coupledwaveguides with the combined on-site quadratic and cubic nonlinearities. Two analyti-cal approximations have been developed, based on the averaging method, which apply,severally, to broad and narrow DESs. Systematic results are obtained by means of numer-ical calculations. Both analytical and numerical methods produce DESs of different types– in particular, with identical or opposite signs of the FF (fundamental-frequency) andSH (second-harmonic) components. The analytical approximation yields very accurateresults, in comparison to their numerically found counterparts, for the most interestingcase of strongly discrete embedded solitons. In that case, the DES family crosses theborder of the spectral propagation band and extends, as a branch of regular solitons, intothe semi-infinite gap. The analytical approximations developed here for the DESs maybe applied to other physically relevant discrete models. Stability of the DESs has beeninvestigated numerically. While it is generally known that DESs are semi-stable objects,which we observed as well here, in the case of FF and SH components having intersitecoupling coefficients with opposite signs, we demonstrate the existence of remarkablystable DESs, which, to the best of our knowledge, has not been reported before.A challenging direction for the development of the present analysis is to perform itfor two-dimensional lattices with the χ (2) : χ (3) nonlinearity, with the aim to seek fortwo-dimensional DESs. Competing interests
The authors declare that they have no known competing financial interests or personalrelationships that could have appeared to influence the work reported in this paper.13
RediT author statementHadi Susanto : Conceptualization, Methodology, Formal analysis, Investigation,Writing - Original Draft, Visualization.
Boris Malomed : Conceptualization, Method-ology, Writing - Review & Editing.
Acknowledgements
HS thanks Mahdhivan Syafwan for fruitful discussions at the early stage of this work.The work of BAM is supported, in part, by Israel Science Foundation through grant No.1286/17. This author appreciates hospitality of Department of Mathematical Sciencesat the University of Essex (UK).
References [1] Y. R. Shen,
The Principles of Nonlinear Optics (Wiley: New York, 1984).[2] Y. N. Karamzin and A. P. Sukhorukov, Nonlinear interaction of diffracted light beams in a mediumwith quadratic nonlinearity – mutual focusing of beams and limitation on efficiency of opticalfrequency converters, Pisma Zh. Exp. Teor. Fiz. , 730 (1974) [JETP Lett. , 338 (1975)].[3] W. E. Torruellas, Z. Wang, D. J. Hagan, E. W. van Stryland, G. I. Stegeman, L. Torner, and C.R. Menyuk, Observation of 2-dimensional spatial solitary waves in quadratic medium, Phys. Rev.Lett. , 5036 (1995).[4] G. I. Stegeman, D. J. Hagan, and L. Torner, χ (2) cascading phenomena and their applications toall-optical signal processing, mode-locking, pulse compression and solitons, Opt. Quantum Electron. , 1691-1740 (1996).[5] C. Etrich, F. Lederer, B. A. Malomed, T. Peschel, and U. Peschel, Optical solitons in media witha quadratic nonlinearity, Prog. Opt. , 483-568 (2000).[6] A. V. Buryak, P. Di Tripani, D. V. Skryabin, and S. Trillo, Optical solitons due to quadraticnonlinearities: from basic physics to futuristic applications Phys. Rep. , 1677-1679 (1995).[8] A. V. Buryak, Y. S. Kivshar, and S. Trillo, Optical solitons supported by competing nonlinearities,Opt. Lett. , 1961-1963 (1995).[9] J. Yang, B. A. Malomed, and D. J. Kaup, Embedded solitons in second-harmonic-generating sys-tems, Phys. Rev. Lett. , 1958-1961 (1999).[10] A. R. Champneys and M. D. Groves, A global investigation of solitary waves solutions to a two-parameter model for water waves, J. Fluid Mech. , 199-229 (1997).[11] J. Fujioka and A. Espinosa, Soliton-like solution of an extended NLS equation existing in resonancewith linear dispersive waves, J. Phys. Soc. Jpn. , 2601-2607 (1997).[12] A. R. Champneys, B. A. Malomed, J. Yang, and D. J. Kaup. “Embedded solitons”: solitary wavesin resonance with the linear spectrum. Physica D , 340-354 (2001).[13] G. Haller, A variational theory of hyperbolic Lagrangian coherent structures, Physica D
0, 574-598 (2011).[14] J. K. Yang, Stable embedded solitons, Phys. Rev. Lett. , 143903 (2003).[15] R. Iwanow, R. Schiek, G. I. Stegeman, T. Pertsch, F. Lederer, Y. Min, and W. Sohler, Observationof discrete quadratic solitons, Phys. Rev. Lett. , 113902 (2004).[16] F. Lederer, G. I. Stegeman, D. N. Christodoulides, D. N. Christodoulides, G. Assanto, M. Segev,and Y. Silberberg, Discrete solitons in optics, Phys. Rep. , 1-126 (2008).[17] K. Yagasaki, A. R. Champneys, and B. A. Malomed, Discrete embedded soliton s , Nonlinearity ,2591-2613 (2005).[18] H. S. Eisenberg Y. Silberberg, R. Morandotti, A. Boyd, and J. S. Aitchison, Diffraction managementPhys. Rev. Lett. , 1863-1866 (2000).[19] T. Pertsch, Y. Silberberg, U. Peschel, and F. Lederer, Anomalous refraction and diffraction indiscrete optical systems, Phys. Rev. Lett. , 093901 (2002).
20] D. Cai, A. R. Bishop, and N. Grønbech-Jensen, Localized states in discrete nonlinear Schr¨odingerequations, Phys. Rev. Lett. , 591-595 (1994).[21] J. H. Dawes and H. Susanto, Variational approximation and the use of collective coordinates, Phys.Rev. E , 063202 (2013).[22] B. A. Malomed and M. I. Weinstein, Soliton dynamics in the discrete nonlinear Schr¨odinger equa-tion, Phys. Lett. A , 91-96 (1996).[23] D. J. Kaup, Variational solutions for the discrete nonlinear Schr¨odinger equation, Math. Comp.Simul. , 322-333 (2005).[24] B. A. Malomed, D. J. Kaup, and R. A. Van Gorder, Unstaggered-staggered solitons in two-component discrete nonlinear Schr¨odinger lattices, Phys. Rev. E , 026604 (2012).[25] D. J. Kaup and B. A. Malomed, Embedded solitons in Lagrangian and semi-Lagrangian System s ,Physica D , 153-161 (2003).[26] M. Syafwan, H. Susanto, S. M. Cox, and B. A. Malomed, Variational approximations for travelingsolitons in a discrete nonlinear Schr¨odinger equation, J. Phys. A: Math. Theor. (7), 075207 (2012).[27] R. Rusin, R. Kusdiantara, and H. Susanto, Variational approximations using Gaussian ansatz, falseinstability, and its remedy in nonlinear Schr¨odinger lattices, J. Phys. A: Math. Theor. (47),475202 (2018).[28] S. Aubry, Breathers in nonlinear lattices: Existence, linear stability and quantization, Physica D , 201-250 (1997)., 201-250 (1997).