Dissipative Nonlinear Josephson Junction of Optical Soliton and Surface Plasmon
DDissipative Nonlinear Josephson Junction of Optical Soliton and Surface Plasmon
Yasa Ek¸sio˘glu, ∗ ¨Ozg¨ur E. M¨ustecaplıo˘glu, and Kaan G¨uven Department of Physics, Ko¸c University, Istanbul 34450 Turkey (Dated: November 7, 2018)We examine the dynamics of a dissipative photonic Josephson junction formed by the weak cou-pling of an optical soliton in a nonlinear dielectric waveguide and a co-propagating surface plasmonalong a parallel metal surface with a linear dielectric spacer. We employ a heuristic model witha coupling function that depends on the soliton amplitude, and consider two phenomenologicaldissipation mechanisms separately: angular velocity dissipation and population imbalance dissipa-tion. In the former dissipation mechanism, the system exhibits phase-slip phenomenon where theodd- π phase modes decay into even- π phase modes. The latter damping mechanism sculptures thephase-space significantly by introducing complex features, among which Hopf type bifurcations arenotable. We show that some of the bifurcation points expand to stable limit cycles for certainregimes of the model parameters. PACS numbers: 42.65.Sf, 42.65.Tg, 03.75.Kk, 03.75.LmKeywords: optical soliton, surface plasmon, josephson junction,nonlinear dynamics, decoherence, nonlinearoptics
I. INTRODUCTION
As the research on the interaction of light with metal-lic structures matures as a well-established technologycalled plasmonics [1–3], the coupling of surface plasmonsto different light sources are being investigated, with themotivation that controlling surface plasmons offer thepotential for developing different types of SP-integratednanophotonic devices [4–6]. In particular, the couplingbetween SPs and confined light modes are widely inves-tigated. A recent proposal is based on the resonant in-teraction between the SPs on a metal surface and co-propagating soliton in a nonlinear dielectric medium [7].Under a classical formulation, this system exhibits richnonlinear dynamical features where the interaction de-pends on the soliton amplitude, as such it may be uti-lized to manipulate the SP propagation. As an addi-tional feature, it has been shown [8] that the system isakin to the bosonic Josephson-junction of Bose-Einsteincondensates [8–12] so that similar and different nonlin-ear Josephson junction features may be realized in thisoptic-plasmonic system.While this was an exciting analogy, surface plasmons aresubject to strong dissipative effects in the host metal andperfect Josephson junction is a too idealized model. Inthe present work, we aim to make a more faithful repre-sentation of the system by taking into account the dis-sipation. Remarkably, we find that the dissipation canbring some benefits if it could be introduced under con-trol. The junction can exhibit stable dynamics in spite ofdissipation. Based on this motivation, we investigate herethe dissipation effects of the coupled soliton-SP systembased on the model introduced in Ref. [7] and then for-mulated as a nonlinear Josephson junction [8]. We note ∗ Electronic address: [email protected] that dissipation effects in the Bosonic Josephson Junction(BJJ) were previously studied, [9] which might be usefulto refer for comparison. We employ a standard dynamicalanalysis of the system in the phase space representation,and investigate how different dissipation mechanisms andother model parameters affect the phase-space landscape.Bifurcations in the system are highlighted which may beinduced in certain regimes of these parameters.
II. MODEL AND THEORETICALFORMULATION
Construction of a classical theoretical model for a cou-pled soliton surface-plasmon system has been addressedrecently in a number of publications. An elegant initia-tive was the heuristic model introduced by Bliokh andcollaborators, [7] in which the interaction between thesoliton and the surface-plasmon is formulated as a cou-pled non-linear/linear oscillator system, with the cou-pling parameter being dependent on the soliton ampli-tude. While being constrained by a number of simplify-ing assumptions, the model predicts stationary coupledmodes of the soliton and the surface-plasmons under fea-sible parameter regimes. In a following work, station-ary and quasi-stationary solutions for the so-called soli-plasmons are investigated under an asymmetric couplingmodel [13]. In this paper, we will employ the heuristiccoupled-oscillator model, albeit with an improved cou-pling parameter adapted from the recent formulations.Our main motivation is to incorporate the correct be-havior for the coupling in the strong soliton and strongsurface-plasmon amplitudes, respectively [13]. We be-gin by recapitulating the heuristic theoretical model inRef. [7], where the interaction between optical solitonsand surface plasmons (SP) is discussed in the context ofan resonant optical system: a soliton propagating par-allel to a metal interface can excite surface-plasmons byits evanescent lateral tail, which, in turn, interacts in a r X i v : . [ phy s i c s . op ti c s ] O c t resonance with the soliton along the propagation. Thesoliton and the SP are assumed to be co-propagating,hence only the spatial dynamics of the propagation isof concern. The nonlinearity of the system is assumedto be confined at a distance, d , from the metal interfacesuch that the surface-plasmon propagation retains linear-ity, and the coupling between the soliton and SP-fieldsis weak. The formulation is based on a 2D coordinatesystem in which y is the propagation direction and x isthe lateral direction.The total electric field of the system is introduced bythe following variational ansatz in which the soliton andSP fields are written in a product form of their respectivelongitudinal ( c p,s ( y )) and transverse amplitudes. E ( x, y ) = c p ( y ) e − κ p x + c s ( y )cosh [ κ s ( x − d )] , (1)Evanescent wavevectors are κ p = (cid:113) k p − k and κ s = k (cid:112) γ/ | c s | . The k is the incoming wave vector whereas k p is the SP wave vector. γ is the nonlinearity parameter ofthe dielectric strip. The longitudinal amplitudes c p,s ( y )obey the coupled spatial propagation wave equations:¨ c p + β p c p = q ( | c s | ) c s , ¨ c s + β s c s = q ( | c s | ) c p . (2)where ¨ c p,s ≡ ∂ c p,s ∂ξ and ξ ≡ ky is the dimensionless prop-agation coordinate. β p = k p /k and β s = 1 + γ | c s | / q ( | c s | ) is the coupling function which we will discuss fur-ther below. The coupled equations are linearized by in-troducing c p,s = C p,s e iξ into Eq. 2 and by employingthe slowly varying amplitude approximation to discardhigher order derivatives. The final set of equations arethen: i ˙ C p = ν p C p − q ( |C s | )2 C p , (3) i ˙ C s = − q ( |C s | ))2 C p + ν s ( |C s | ) C s . (4) ν p ≡ β p − (cid:28) ν s ≡ β s − (cid:28) q ( |C s | ) (cid:39) |C s | e − σkd √ γ |C s | (5)Here, σ is an overall constant which we introduce as aparameter to account for the other constants (e.g. per-mittivities, propagation constants). The key parametersare preserved in the functional form. In the strong solitonamplitudes ( |C s | ≈ |C s | since the soliton lat-eral profile in Eq. 1 would be wider resulting in a largeroverlap.The analogy between this model and the bosonicJosephson junction (BJJ) dynamics [10, 12] is revealedwhen we further substitute C s,p = C s,p e iφ s,p into Eq. 4and introduce a new variable set consisting of the frac-tional population imbalance Z = ( | C s | − | C p | ) /N andthe relative phase difference between the soliton and theSP φ = φ s − φ p :˙ Z = − q ( Z ) (cid:112) − Z sin φ, (6)˙ φ = Λ Z + ∆ E + q ( Z ) Z √ − Z cos φ, (7)where Λ = γN , is the nonlinearity parameter, and ∆ E ≡ Λ − ν p parametrizes the asymmetry between the solitonand SP states (similar to the asymmetry between thewells of a double well system). N = ( | C s | + | C p | ) isconstant for the isolated system. The coupling function q takes the form q ( Z ) (cid:39) (cid:114) (1 + Z )2 e − σkd √ Z ) . (8)We stress that the Z dependence makes an inherently dy-namic coupling as opposed to constant (or externally tun-able) coupling parameter present in BJJ-systems. Withthis formulation, the coupled soliton SP system can bedescribed as a nonlinear Josephson junction.The range of the model parameters are chosen to bearound the resonant coupling regime [7]. In the calcula-tions presented below, we fix the dimensionless separa-tion parameter kd = 6, and choose the other parametersbetween Λ (cid:39) . − .
03, ∆ E (cid:39) -0.4 − kd values would confine the inter-action to smaller soliton amplitudes and eventually resultin the decoupling of the soliton and the SP completely.Small kd values would not be applicable in the weak cou-pling approximation. III. ANGULAR-VELOCITY DEPENDENTDISSIPATION
Owing to the analogy to the BJJ models [9], the firstdissipation mechanism we consider is about the incoher-
FIG. 1: (Colour online)(a)The Z − φ phase-space of the dissipationless SP-soliton coupled system with kd = 6, Λ = 0 .
01, ∆ E = − . η = 0 .
2. Phase slip occurs between the φ = − π and φ = 0 modes. ent exchange of photons between the soliton and the SP.This is introduced by the term η ˙ φ in the following equa-tions, ˙ Z = − q ( Z ) (cid:112) − Z sin φ − η ˙ φ, (9)˙ φ = Λ Z + ∆ E + q ( Z ) Z √ − Z cos φ, (10)Under coordinate reversal ξ → − ξ , we have Z → Z and φ → − φ , hence, the dissipative term η ˙ φ appearing inthe equation for ˙ Z represents an irreversable process cor-rectly. The phase-space of the dissipationless and dissipa-tive system are shown in Fig.1(a) and 1(b), respectively.The gradient distribution is given with color-coded mag-nitudes and typical trajectories are indicated by whitecurves. The phase space is characterized by the zero-(even) phase and π -(odd) phase modes, where the criti-cal points are located. There is always one stable pointat φ = 0. At φ = π , one (stable) to three (two stable,one saddle) critical points may emerge depending on thevalues of the model parameters Λ , kd , and ∆ E [8]. InFigure 1(a), sample trajectories (white curves) show an-harmonic closed orbits around stable critical points andopen trajectories. The color indicates that the gradientsget steeper in the SP-dominant population imbalance re-gion ( Z < Z =0 , ˙ φ = 0 condition, the ˙ φ -dependent dissipation does notalter the location of the critical points in the phase spacebut the eigenvalues of the Jacobian acquire nonzero realparts such that the stable critical points at φ = π be-comes source and the critical point at φ = 0 becomes asink. While the trajectories flow towards the sink points,the phase-slip phenomenon is observed from φ =odd- π -modes to φ =even- π modes, as reported for the BJJsystems [9]. Depending on the value of the asymme-try parameter, ∆ E , the phase-slip occurs with a nonzerochange in the averaged population imbalance < Z > . With a larger (smaller) dissipation constant, the typ-ical phase-diagram depicted in Fig.1(b) gets compressed(extended) along the phase axis but the features de-scribed above remain the same. Thus, while the tran-sient behavior depends on η , the final state of the systemis independent: φ s = φ p and a fixed value of Z . IV. POPULATION IMBALANCE DISSIPATION
In the following, we introduce the set of equations thatincorporate a dissipation proportional to the populationimbalance, Z . In the context of BJJ systems, such a dis-sipation arises from the finite lifetime of excited states oftwo condensates in different hyperfine levels in a singleharmonic trap connected by tunnelling transitions [14].In the analogous mechanical system of momentum short-ened pendulum, this corresponds to a damping propor-tional to the angular momentum variable [15], which weconsider to be applicable to the soliton-SP system phe-nomenologically.˙ Z = − q ( Z ) (cid:112) − Z sin φ − ζZ, (11)˙ φ = Λ Z + ∆ E + q ( Z ) Z √ − Z cos φ. (12)The first thing we note here is that this dissipationterm can modify the critical points in the phase-spacesignificantly: their existence, location and correspondingJacobian eigenvalues. For determining the critical points,using ˙ Z = 0 , ˙ φ = 0 and the identity sin ( φ ) + cos ( φ ) = 1we obtain the following root equation in Z , Z ζ +(1 − Z )[(∆ E +Λ Z ) (1 − Z ) − q ( Z ) Z ] = 0 (13)Figure 2 is a plot this equation for fixed kd = 6 , Λ = 0 . E = 0, and (b) ∆ E = − . ζ . Blue FIG. 2: (Colour online) (a)Plot of Eq. 13 for kd = 6 , Λ = 0 .
01 and ∆ E = 0 at various values of the dissipation constant ζ . Roots of f ( Z ) are thecritical points in the phase-space,(b)Plot of Eq. 13 for kd = 6 , Λ = 0 .
01 and ∆ E = − . ζ . Rootsof f ( Z ) are the critical points in the phase-space. curve is plotted for dissipationless case ( ζ = 0) for refer-ence. For ∆ E = 0, Z = 0 is always a critical point as allcurves in Fig. 2(a) are tangent to Z = 0 at this point,irrespective of dissipation. This can also be inferred fromEq. 13. Critical points close to Z ≈ ζ ≈ − ). For ∆ E = − . ζ = 0 . ζ > .
01. Table I classifies the critical pointsat different dissipation regimes for ∆ E = − . TABLE I: Classification of fixed points at different dissipationregimes for Fig. 2(b) ζ I II III IV V VI0 <ζ< . × − saddle sso ssi saddle ssi saddle8 . × − ≤ ζ< .
016 saddle sso ssi saddle - -0 . ≤ ζ< .
027 saddle sso - - - -0 . ≤ ζ - - - - - -sso: spiral source, ssi:spiral sink According to the typical ranges of the dissipation pa-rameter we now investigate the behavior of the systemin the phase-space. Figure 3(a) is plotted for kd = 6,Λ = 0 .
01, ∆ E = 0, and ζ = 10 − . These values corre-sponds to the cyan curve in Fig. 2(b). Z = 0 is a spiralsink at φ = 0 and a saddle point at φ = π . Another spi-ral sink is located at Z = − . , φ = 0 . π . The saddlepoint at Z = − . φ = 0 . π is not visible at this scale.In Figure 3(b), same parameters as in Fig. 3(a) areused but ∆ E = − . ζ = 0 .
02 for whicha larger limit cycle is observable. Finally, ζ ≈ .
028 isthe onset of the regime in which all trajectories flow intoa unique stable cycle and there are no critical points inthe phase space. This behavior is similar to that of aJosephson junction between two superconductors, withbias current above the critical current, I ≡ I B I C > φ and an-gular velocity ˙ φ ≡ y satisfy the dimensionless equations˙ φ = y, (14)˙ y = I − sin( φ ) − αy. (15)where α is the damping parameter. Note that for I >
A. Codimension-1 Bifurcations
In order to investigate the bifurcation dynamics, wefirst consider codimension one bifurcations which can beinduced by varying a single parameter of the model. Weemploy a numerical continuation software called Mat-cont [17], which computes how a critical point advances inthe phase space by varying model parameters, and char-acterizes the point along the continuation according toits Jacobian. We begin by the continuation with respectto the dissipation constant ζ . The fixed parameters are kd = 6, Λ = 0 .
01, and ∆ E = − . FIG. 3: (Colour online) Phase-space plot of coupled soliton-SP system with population imbalance dissipation. Model parameters are kd = 6 , Λ =0 . ζ = 0 .
01. (a) ∆ E = 0, (b) ∆ E = − . FIG. 4: (Colour online)Phase-space plot of coupled soliton-SP system with population imbalance dissipation. Model parameters are kd = 6 , Λ =0 .
01, ∆ E = − . ζ = 0 .
02. Exterior and interior trajectories shown for left and right limit cycles, respectively. (b) ζ = 0 . of limit cycles in Z axis, emanating from the Hopf-point.Figure 5 (b) shows the limit cycles in the Z − φ − ζ plane.Interestingly, the limit cycles get bigger in the φ − Z planefor stronger dissipation, before they are terminated by alimit point cycle. For reference, the limit cycles in Fig-ures 3 and 4 are in the range of Fig. 5(b). TABLE II: Bifurcation points shown in Fig. 5(a)Label ∆
E φ ∗ Z ∗ a H 0 . − . − . − . . − . − . − . . − . − .
495 -LP: Limit Point (fold),H:Andronov-Hopf, NS:NeutralSaddle, a is the normal form coefficient for LP and firstLyapunov coefficient for H.
Next, the continuation is performed in the asymmetryparameter ∆ E . Figures 6(a) and 6(b) show these resultsfor ζ = 0 . E >
E > Z − ∆ E planefor different values of Λ in Fig. 7 and for ζ in Fig. 8. Equi-librium points are periodic in the Z − φ plane, thus, eachcontinuation curve closes onto itself in these parametricplanes. The continuation curves extend in a finite regionof ∆ E − Z plane, where the range also depends on othermodel parameters.In Figure 7, we depict the effect of nonlinearity (Λ) onthe continuation curve. In view of Eq. 8 we note that,for a given soliton amplitude, increasing nonlinearity con-fines the soliton more tight around its propagation axisand decreases the coupling between the soliton and theSP. Thus, while the continuation curve extends in thewhole Z axis for small Λ (blue curve) it shrinks and gets FIG. 5: (Colour online)(a) Continuation of the stable point as a function of dissipation constant. Model parameters are kd = 6 , Λ = 0 . E = − . Z -axis. A limit point cycle terminates the limit cycles at ζ ≈ .
01, (b) Limit cycles plotted in the Z − φ − ζ parametricphase-space. FIG. 6: (Colour online) (a)Continuation of the stable point as a function of ∆ E . Model parameters are kd = 6 , Λ = 0 . ζ = 0 . Z − φ − ∆ E parametric phase-space. TABLE III: Bifurcation points shown in Fig. 6(a)Label ∆
E φ ∗ Z ∗ a LP − .
340 0 . − .
997 3 . − .
148 0 . − . − . − . . − . − . . − . − .
139 0 . . − .
311 0 . − . . . − . − . .
167 0 . − . − . .
360 0 . − . − . confined to the range close to Z = − E value when | ∆ E | is small (e.g. < . .
003 only LP occur. For Λ = 0 .
03 LP and H bifurca-tions are present but the first Lyapunov exponent is posi-tive (i.e. subcritical Hopf bifurcations). For the interme-diate values between these extremes, we observe that thecontinuation curve hosts multiple H-points which appearin pairs in negative and positive sides of ∆ E axis. Theseare supercritical, and as it was shown in Fig. 6 stablelimit cycles are spawn between them.When the dissipation strength ζ increases, the contin-uation curve in the ∆ E − Z plane shrinks as shown inFig. 8. Limit point and Hopf bifurcations decorate thecontinuation curves. The Hopf points close to Z = − ζ > .
028 sinceno critical point is present.We recall that ∆ E = Λ − ν p is the difference be-tween the (soliton-characteristic) nonlinearity and the(SP-characteristic) propagation constant. While not ob-vious in these terms, it is analogous to the asymmetry FIG. 7: (Colour online) ∆ E continuation of the stable point for dif-ferent values of Λ, with kd = 6, ζ = 0 . FIG. 8: (Colour online) ∆ E continuation of the stable point for dif-ferent values of ζ , with kd = 6, Λ = 0 .
01. Hopf (H) and limit pointbifurcations are labeled on the continuation curves. Inset shows en-largement around the ∆ E = 0 , Z = 0 region. Blue curve correspondsto the yellow curve in Fig. 7 between the wells of the double well potential in BJJsystems, and thus, sizes the asymmetry between the soli-ton state and the SP state of the Josephson junction. InFigures 7 and 8, we observe that for large magnitudes of∆ E (i.e. irrespective of its sign), the equilibrium pointslie in the SP-dominant ( − < Z <<
0) region. Solitondominant equilibrium points occur under weak nonlin-earity and weak dissipative conditions only.
B. Codimension-2 Bifurcations
Codimension-2 bifurcations require the simultaneousvariation of two parameters. The blue curve in Fig. 9shows the continuation an equilibrium point in the ζ − ∆ E parameter space. We observe three Bogdanov-Takens(BT) bifurcation points and a cusp point. BT-bifurcationimplies the existence of three codimension-one bifur-cations nearby, namely, a saddle-node bifurcation, anAndronov-Hopf bifurcation and a homoclinic bifurcation.From the BT points, the continuation curves of the Hopfpoint and limit point can be obtained respectively. The orange curve is the Hopf curve passing through the BTpoint close to the cusp point. The green curve is theHopf curve passing through the lower BT point. EachHopf curve has two generalized Hopf-points when theyintersect the ζ = 0 axis. The black curve is the continu-ation of the limit point from the upper BT point on thegreen Hopf curve.This codimension-2 diagram indicates the regions of ζ and ∆ E in which stable limit cycles can be generated inthe system. In Figure 6 it was shown that for fixed ζ stable limit cycles are spawn between two Hopf points.Here, this corresponds to the region between the orangeand the green curves. The Cusp point on the blue curvemarks the value of ζ = 0 .
028 beyond which there areno critical points for Λ = 0 . E = − . ζ − ∆ E plane and φ − Z plane and theirnormal form coefficients are listed. BT points have twonormal form coefficients with respect to two varying pa-rameters. Cusp point is characterized by a single normalform coefficient. FIG. 9: (Colour online) Codimension-2 bifurcation diagram in the ζ − ∆ E plane. Model parameters are kd = 6 , Λ = 0 .
01. Blue curveshows the continuation of an equilibrium point. Orange and greencurves are the Hopf continuation curves from corresponding BT points.Black curve is the limit point continuation curve. BT: Bogdanov-Takens point, C: Cusp point. GH: Generalized Hopf point. Inset showsthe zoomed view around the ∆ E = 0 axis. V. CONCLUSION
In this work, we investigated the dynamical featuresof a dissipative nonlinear Josephson junction which isformed by a paraxial optical soliton surface-plasmon cou-pled system. The most notable part of the heuristicmodel implemented here is the nonlinear nature of thecoupling mechanism. For a dissipation proportional tothe angular velocity, the system decays into stable fixedpoints with constant relative phase and constant averagepopulation imbalance. A particular feature observed hereis the phase-slip phenomenon, where the odd- π modes ofthe dissipationless system decay into even- π modes. Fordissipation proportional to the population imbalance, su-percritical Andronov-Hopf bifurcation can occur which isnotable in that stable oscillatory modes between the soli- TABLE IV: Codimension-2 bifurcations in the ζ − ∆ E parameter space (blue curve in Fig 9)Label ζ ∆ E φ ∗ Z ∗ ( a, b ) (1) , c (2) BT 0 3 . × − − . − .
139 (1 . × − , − . × − )BT 0 .
026 1 . × − − . − .
434 (1 . × − , − . × − )CP 0 .
028 3 . × − − . − . − . . − . − . − .
946 ( − . , . (1) normal form coefficients of BT. (2) normal form coefficient of CP. ton and surface-plasmon fields can be maintained againststrong dissipation. The heavily damped limit resemblesthe behavior of Superconducting Josephson junctions inwhich the system decays rapidly into a unique oscillat-ing state. In view of this analysis, we propose that thenonlinear Josephson junction between the optical solitonand surface plasmons can exhibit interesting dynamicalfeatures which may be exploited further based on morerigorous theoretical models. Acknowledgments
We acknowledge the support by the Science and Tech-nology Research Council of Turkey (TUBITAK) underthe project no: 111T285. K. G¨uven acknowledges thesupport by the Turkish Academy of Sciences. [1] A. V. Zayats, I. I. Smolyaninov and A. A. Maradudin,Phys.Rep. , 131 (2005).[2] S. A. Maier
Plasmonics: fundamentals and applications (Springer, New York, 2007).[3] A. R. Davoyan, I. V. Shadrivov and Y. S. Kivshar, OpticsExpress , 21732 (2009).[4] (Ed.s) M. L. Brongersma M L and P. Kik Surface Plas-mon Nanophotonics (Springer, Netherlands, 2007).[5] R. H. Ritchie, Phys. Rev. , 874 (1957).[6] K. Y. Bliokh, Y. P. Bliokh, V. Freilikher, S. Savelev andF. Fori, Rev. Mod. Phys. , 1201 (2008).[7] K. Y. Bliokh, Y. P. Bliokh and A. Ferrando, Phys. Rev.A , 041803R (2009).[8] Y. Ek¸sio˘glu, . E. M¨ustecaplıo˘glu, K. G¨uven, Phys. Rev.A , 033805 (2011).[9] I. Mariono, S. Raghavan, S. Fantoni, S. R. Shenoy andA. Smerzi, Phys. Rev. A , 487 (1999).[10] A. Smerzi, S. Fantoni, S. Giovanazzi and S. R. Shenoy,Phys. Rev. Lett. , 4950 (1997).[11] S. Giovanazzi, A. Smerzi and S. Fantoni, Phys. Rev. Lett. , 4521 (2000).[12] S. Raghavan, A. Smerzi, S. Fantoni and S. R. Shenoy, Phys. Rev. A , 620 (1999).[13] C. Milian, D. E. Ceballos-Herrera, D. V. Skyrabin andA. Ferrando, arXiv:1205.3182v1, (2012).[14] C. Lee, W. Hai, L. Shi, X. Zhu, K. Gao, Phys. Rev. A , R31 (1999).[16] Steven H. Strogatz, Nonlinear Dynamics and Chaos withApplications to Physics, Biology, Chemistry and Engi-neering (Perseus Books Publishing L.L.C, 1994).[17] MatCont is a Matlab software project for the numer-ical continuation and bifurcation study of continuousand discrete parameterized dynamical systems devel-oped under the supervision of W. Govaerts and Yu.A.Kuznetsov(http://sourceforge.net/projects/matcont/).[18] J.Hale and H.Kocak,
Dynamics and Bifurcations (Springer-Verlag, 1991).[19] D.K Arrowsmith and C.M. Place,