ABSTRACT

In this paper, a methodology for the numerical location of a global point-to-point (P2P for short) homoclinic asymptotically connecting orbit is applied to a modified version of Shimizu-Morioka system, which models a semiconductor laser. This type of global bifurcation can be considered as a stylized mathematical description of self-pulsation in this laser type, associated with saturation. The location is achieved by use of a custom algorithm based on the method of orthogonal collocation on finite elements with fourth order boundary conditions, constructed through scale order approximations. The effectiveness of the algorithm and the superiority of high-order boundary conditions over the widely used first order ones are justified throughout the obtained graphical results.

Keywords:

Homoclinic Connections, Orthogonal Collocation on Finite Elements, Numerical Continuation, Limit Cycles

1. Introduction

The numerical location of global asymptotic orbits is often proved to be a computationally demanding task, even in low dimensional systems. However, the recent improvement of hardware capabilities and symbolic mathematical software allows researchers to locate global asymptotic orbits numerically more easily. Homoclinic point-to-point (P2P for short) connecting orbits arise in various applications where hysteresis and saturation phenomena are present. These orbits are considered to be separatrices in the nonlinear state space in 2D conservative ODEs (ordinary differential equations), since they divide the phase space into qualitatively different regions of motion; one region of periodic solutions and one of non-periodic ones, respectively. In that sense, homoclinic orbits can be perceived as the limit of periodic solutions, which is a periodic orbit, the fundamental period of which tends to infinity, while the orbit itself remains bounded. From a different perspective, a homoclinic P2P orbit can be regarded as the result of the collision of a limit cycle and a fixed point.

In the present paper, an algorithm for the numerical computation of global homoclinic P2P orbits is applied to a three-dimensional differential dynamical system, be it a modified Shimizu-Morioka system [1] , [2] modelling a semiconductor laser. The aforementioned algorithm has been carried out by means of the well-known method of orthogonal collocation on finite elements (widely used in the famous software AUTO86 (see [3] ) and MATCONT (see [4] , [5] )) and it has already been successfully applied by the authors to a Lorenz system, as well as to a three-species food chain model with group defence ecosystem (see [6] ). In the version applied herein, Lagrange polynomials have been used as the basis polynomials, instead of the Legendre polynomials implemented in [6] (see Section 2, Equation (4)). Hence, the collocation points associated with each time subinterval (see [6] , Figure 1) are placed at the Gauss points, which are the roots of the maximum degree Legendre polynomial relative to this subinterval.

The respective procedure involves the evaluation of high order boundary conditions (BC from now on). Note that both the well-known and widely used techniques of projection BC and the method of eigenvectors (see [7] ) are approximations of the first order. Thus the application of those techniques presupposes good initial data for the successful computation of the orbits of interest and this is becoming harder and harder to achieve as the number of state variables together with the number of active parameters increase. Moreover, it is well known that the projection method converges exponentially with the increase of

the truncation interval, which in turn, however, can increase the computational time (mainly in ordinary PCs with low to moderate CPU power such as an Intel Core i7 870). So, the use of high order BC can be proved useful in cases like that. An appropriate combination of the multiple scales approximation method with that of successive approximations leads to the construction of a technique for the determination of high order BC. The idea for this combination comes from Deprit and Henrard [8] , Bennett [9] and the relevant references therein. Also, Hassard [10] presented the idea to use high order BC instead of the projection ones.

The theoretical background of the computational procedure is given in [6] (Sections 2 and 3). More precisely, in the paper mentioned above, the defining equations and important formulae associated with the algorithm implemented are presented throughout Section 2, while a brief description of the procedure concerning the derivation of high order BC is presented in Section 3. In the same Section, the determination of the number of control parameters is treated. This concerns the number of the necessary active parameters of the system for which the homoclinic orbit is an isolated, structurally stable phenomenon. By the analysis developed in [6] (Section 3), for the system presented herein this number has also been calculated equal to one, so we have chosen one active parameter, a (see Section 2 below). Thus, in Section 2 of the present paper, the algorithm is applied to a modified version of Shimizu-Morioka system, modelling a semiconductor laser. After a brief review of the theory associated with the type of lasers under consideration, where the significance of homoclinic orbits is mentioned, as well, an equilibrium analysis is performed and the scale order approximations (up to the fourth order) of the system together with the respective outgoing solution are obtained. Both the latter and the similarly constructed incoming solution lead to the extraction of high (fourth) order BC, the use of which together with the method of orthogonal collocation on finite elements results in the homoclinic orbit under consideration. The analysis is carried out with the aid of Mathworks Matlab and the symbolic engine of Maplesoft Maple, concluding with the associated graphical results.

2. Application to a Modified Version of Shimizu-Morioka System

The method of orthogonal collocation on finite elements with fourth order BC has been applied to a modified version of Shimizu-Morioka system [1] , [2] , for which the location of a homoclinic connection at the fixed point $E_0\left(0,0,0\right)$ has been carried out. The system is described by the following differential equations

$\begin{array}{l}\stackrel{\dot{}}{x}=y\\ \stackrel{\dot{}}{y}=x-qy-xz\\ \stackrel{\dot{}}{z}=-bz+b{x}^2\end{array}$ (1)

and the connection between its parameters and the Lorenz ones is presented in [11] . This version of the system deals with the modelling of a homogeneously broadened single-mode laser (see [12] and the relevant references therein). By scaling $t=\tau /b,y=bY$ , setting $q=\left(1-b^2\right)/b$ and then $b^2=1/\left(1+a\right)$ (see [2] ), (1) becomes

$\begin{array}{l}\stackrel{\dot{}}{x}=y\\ \stackrel{\dot{}}{y}=\left(a+1\right)x-ay-\left(a+1\right)xz\\ \stackrel{\dot{}}{z}=-z+x^2\end{array}$ (2)

where dot denotes differentiation with respect to $\tau$ and we have resubstituted $\tau$ with t and Y with y. The state variables $x,y,z$ describe the small amplitude dynamics of the Laser with Saturable Absorber (LSA) [13] , and more precisely the amplitude of electric field intensity E, the amplitude of atomic polarization P and the atomic polarization differences in absence of the laser field, respectively. Moreover, a is the active (control or bifurcation) parameter.

2.1. Brief Theory and Significance of Homoclinic Orbits

The LSA lasers exhibit sustained laser oscillations consisting of pulse trains of really short high intensity and high frequency laser output (also known as passive Q-switching or self-pulsing behaviour). Such behaviour has been both theoretically and experimentally confirmed for CO₂ lasers, microchip solid state lasers etc. Semiconductor lasers exhibit rates of high repetition ranging from hundreds of MHz to almost 1/10 GHz, while they have useful applications in telecommunication and in optical data storage using CD and DVD systems. Some other applications include the optical feedback noise reduction in semiconductor injection lasers and optical timing extraction by injection locking of self-pulsing optical oscillators. The phenomenon of self-pulsation is a result of the nonlinear interaction of the slowly responding amplifying and absorbing media and the fast response of the electric field in lasers driven by a constant pumping power. Actually, the basic mechanism responsible for the generation of these oscillations starts as the laser is turned on and the amplifying medium, the gain, is excited to a sufficiently high level via some type of pumping process. The absorber absorbs the free photons in the laser and thus the intensity of the electric field remains low and the saturation of the gain goes on. As soon as the absorbing medium saturates, the usual laser process starts with a strongly excited gain causing a high electric field intensity and thus truly enhanced output power. During this process both the gain and the absorber return to ground state and the process starts all over again. In semiconductor lasers, this can produce a pulse train with a typical frequency of the order of several GHz. The homoclinic point-to-point connecting orbit arising in the system of interest, can be considered as a mathematical representation of the aforementioned high power self-pulsation. A characteristic of lasers with obvious practical importance is that by changing either the material or the pump power, the qualitative behaviour of the laser beam can change dramatically, that is a large variety of local and global bifurcations can occur.

However, the presence of quite different time scales makes numerical simulations of lasers both challenging and time consuming, so the use of appropriate time scaling transformations can be useful, sometimes.

2.2. Computation of Equilibria

By means of (2) we easily find the three equilibria of the system, be them $E_0\left(0,0,0\right),E_{\pm }\left(\pm 1,0,1\right)$ . Also, the Equation (2) are invariant under the transformation $x\to -x,y\to -y$ , so its orbits are symmetric with respect to z-axis and we can restrict our analysis to $x>0$ . Moreover, the Jacobian matrix associated with (2), is

$J=\left(\begin{array}{ccc}0& 1& 0\\ a+1-\left(a+1\right)z_0& -a& -\left(a+1\right)x_0\\ 2x_0& 0& -1\end{array}\right)$ (3)

where $\left(x_0,y_0,z_0\right)$ represents the equilibrium under consideration (the nonlinearity condition for higher order terms is obviously valid). Thus, for $E_0$ we have

${J|}_{E_0}=\left(\begin{array}{ccc}0& 1& 0\\ a+1& -a& 0\\ 0& 0& -1\end{array}\right)$ (4)

giving the eigenvalues ${\lambda }_1=-1,{\lambda }_2=1,{\lambda }_3=-a-1$ , hence $E_0$ is a saddle. Furthermore, with regard to $E_{+}$ the Jacobian matrix becomes

${J|}_{E_{+}}=\left(\begin{array}{ccc}0& 1& 0\\ 0& -a& -a-1\\ 2& 0& -1\end{array}\right)$ (5)

with characteristic equation

${\lambda }^3+\left(a+1\right){\lambda }^2+a\lambda +2\left(a+1\right)=0$ (6)

By setting $p_0=2\left(a+1\right),p_1=a,p_2=a+1$ , then according to Routh-Hurwitz criterion, as long as the following relations:

$\begin{array}{l}{p}_0>0\Rightarrow a+1>0\\ p_1>0\Rightarrow a>0\\ D_2=p_1{p}_2-p_0>0\Rightarrow \left(a-2\right)\left(a+1\right)>0\end{array}$ (7)

hold together, that is for $\alpha >2$ (only positive values of the parameter have a physical meaning for our model), the equilibrium $E_{+}\left(1,0,1\right)$ is locally asymptotically stable, since in that case the eigenvalues have negative real part. Also, according to Liu criterion [14] , in the case where the first two inequalities of (7) are valid, then the following equality:

$p_1{p}_2-p_0=0\Rightarrow \left(a-2\right)\left(a+1\right)=0\Rightarrow a=a_0=2$ (8)

determines a Hopf bifurcation, if in addition $\text{d}{D}_2\left(a_0\right)/\text{d}a\ne 0$ also holds, which is equivalent to the transversality condition. Since the critical first Lyapunov coefficient is evaluated $l_1\left(a_0\right)\simeq -5.183052\times {10}^{-2}$ (by means of the corresponding normal form formulae included in the algorithm presented in [15] ), a nondegenerate (codimension 1) supercritical Hopf bifurcation occurs, that is stable limit cycles bifurcate from the unstable equilibrium ( $a<2$ ). Then, the strategy of analysis is to start from a limit cycle, which bifurcates from the equilibrium $E_{+}$ , and numerically continue it with respect to its period and the necessary (active) parameter towards the direction of the increase of the period, until it gets close enough to the saddle fixed point of interest, $E_0$ , associated with the homoclinic orbit, and also the active parameter remains practically unchanged. So, as the largest cycle can be considered as a good initial approximation of the homoclinic orbit, the main computation of the homoclinic orbit can be carried out. The numerical continuation herein has been carried out by means of a custom algorithm of sequential numerical continuation based on the method of orthogonal collocation on finite elements and the corresponding integral phase condition (see [6] , Section 2, Equation (8)). The corresponding numerically continued cycles are presented in Figure 1.

2.3. Application of High Order Boundary Conditions

Let us now apply the aforementioned high order BC. Assuming the solutions of the differential system under consideration can be approximated by

$x\left(t\right)\approx {\displaystyle \underset{i=1}{\overset{k}{\sum }}{\epsilon }^i{x}_i}\left(t\right),y\left(t\right)\approx {\displaystyle \underset{i=1}{\overset{k}{\sum }}{\epsilon }^i{y}_i}\left(t\right),z\left(t\right)\approx {\displaystyle \underset{i=1}{\overset{k}{\sum }}{\epsilon }^i{z}_i}\left(t\right)$ (9)

where $\epsilon$ is the orbital parameter and k is the maximal order of approximation. Then, by substituting the expressions of (9) into (2) (the Taylor expansion with respect to $E_0$ coincides with the system itself) and equating the same order terms, we obtain the respective linear (with respect to the variables corresponding to the j-order, $j=1,\cdots ,k$ ) systems. In terms of the present analysis the maximal order of approximation has been chosen to be $k=4$ . Then, each system is solved and the solution is substituted in the subsequent, higher order systems (method of successive approximations). Let us present these systems.

1st order of approximation

$\begin{array}{l}{\stackrel{\dot{}}{x}}_1=y_1\\ {\stackrel{\dot{}}{y}}_1=\left(a+1\right)x_1-a{y}_1\\ {\stackrel{\dot{}}{z}}_1=-z_1\end{array}$ (10)

2nd order of approximation

$\begin{array}{l}{\stackrel{\dot{}}{x}}_2=y_2\\ {\stackrel{\dot{}}{y}}_2=\left(a+1\right)x_2-a{y}_2-\left(a+1\right)x_1{z}_1\\ {\stackrel{\dot{}}{z}}_2=-z_2+x_1^2\end{array}$ (11)

3rd order of approximation

$\begin{array}{l}{\stackrel{\dot{}}{x}}_3=y_3\\ {\stackrel{\dot{}}{y}}_3=\left(a+1\right)x_3-a{y}_3-\left(a+1\right)\left(x_1{z}_2+x_2{z}_1\right)\\ {\stackrel{\dot{}}{z}}_3=-z_3+2x_1{x}_2\end{array}$ (12)

4th order of approximation

$\begin{array}{l}{\stackrel{\dot{}}{x}}_4=y_4\\ {\stackrel{\dot{}}{y}}_4=\left(a+1\right)x_4-a{y}_4-\left(a+1\right)\left(x_1{z}_3+x_2{z}_2+x_3{z}_1\right)\\ {\stackrel{\dot{}}{z}}_3=-z_3+2x_1{x}_3+x_2^2\end{array}$ (13)

By means of the procedure described in [6] (see Section 3, p. 558), we arrive at the approximations of both the outgoing (locally asymptotically unstable) vector solution and the incoming (locally asymptotically stable) one. These approximations can be extracted with the aid of a symbolic computational package, such as Mathematica or Maplesoft Maple (which offers direct integration with Mathworks Matlab), as the calculations can be lengthy even for low dimensional systems. We present below the solutions associated with the outgoing vectors:

1st order approximation

$\begin{array}{l}{x}_1\left(t\right)=c_1{\text{e}}^t+c_2{\text{e}}^{-\left(a+1\right)t}\\ y_1\left(t\right)=c_1{\text{e}}^t+c_2\left(-a-1\right){\text{e}}^{-\left(\alpha +1\right)t}\\ z_1\left(t\right)=c_3{\text{e}}^{-t}\end{array}$ (14)

where $c_1,c_2,c_3$ denote the integration constants (from now on $c_i,i=1,2,3,4,\cdots$ will denote integration constants unless stated otherwise). By setting $c_2=c_3=0$ we get the first order approximation of the outgoing solution vector

$\begin{array}{l}{x}_{1out}=c_1{\text{e}}^t\\ y_{1out}=c_1{\text{e}}^t\\ z_{1out}=0\end{array}$ (15)

2nd order approximation

$\begin{array}{l}{x}_2\left(t\right)=c_4{\text{e}}^t+c_5{\text{e}}^{-\left(a+1\right)t}\\ y_2\left(t\right)=c_4{\text{e}}^t-c_5\left(a+1\right){\text{e}}^{-\left(\alpha +1\right)t}\\ z_2\left(t\right)=c_6{\text{e}}^{-t}+\frac{1}{3}{c}_1^2{\text{e}}^t\end{array}$ (16)

and setting $c_4=c_5=c_6=0$ we take the second order approximation

$\begin{array}{l}{x}_{2out}\left(t\right)=0\\ y_{2out}\left(t\right)=0\\ z_{2out}\left(t\right)=\frac{1}{3}{c}_1^2{\text{e}}^{-t}{\text{e}}^{2t}\end{array}$ (17)

3rd order approximation

$\begin{array}{l}{x}_3\left(t\right)=c_7{\text{e}}^t+c_8{\text{e}}^{-\left(a+1\right)t}-\frac{a+1}{6\left(a+4\right)}{c}_1^3{\text{e}}^{3t}\\ y_3\left(t\right)=c_7{\text{e}}^t-c_8\left(a+1\right){\text{e}}^{-\left(a+1\right)t}-\frac{a+1}{2\left(a+4\right)}{c}_1^3{\text{e}}^{3t}\\ z_3\left(t\right)=c_9{\text{e}}^{-t}\end{array}$ (18)

The third order approximation is obtained by setting $c_7=c_8=c_9=0$ ,

$\begin{array}{l}{x}_{3out}=-\frac{a+1}{6\left(a+4\right)}{c}_1^3{\text{e}}^{3t}\\ y_{3out}=-\frac{a+1}{2\left(a+4\right)}{c}_1^3{\text{e}}^{3t}\\ z_{3out}=0\end{array}$ (19)

4th order approximation

$\begin{array}{l}{x}_4\left(t\right)=c_{10}{\text{e}}^t+c_{11}{\text{e}}^{-\left(a+1\right)t}\\ y_4\left(t\right)=c_{10}{\text{e}}^t-\left(a+1\right)c_{11}{\text{e}}^{-\left(a+1\right)t}\\ z_4\left(t\right)=c_{12}{\text{e}}^{-t}-\frac{a+1}{15\left(a+4\right)}{c}_1^4{\text{e}}^{4t}\end{array}$ (20)

and setting $c_{10}=c_{11}=c_{12}=0$ we get the fourth order approximation of the outgoing solution

$\begin{array}{l}{x}_{4out}\left(t\right)=0\\ y_{4out}\left(t\right)=0\\ z_{4out}\left(t\right)=-\frac{a+1}{15\left(a+4\right)}{c}_1^4{\text{e}}^{4t}\end{array}$ (21)

So, by substituting the aforementioned formulae of $\left(x_j^{out},y_j^{out},z_j^{out}\right),j=1,2,3,4$ for $\left(x_j,y_j,z_j\right)$ in (9), we arrive at the outgoing solution up to the fourth order:

$\begin{array}{l}{x}_{out}=\epsilon c_1{\text{e}}^t-\frac{a+1}{6\left(a+4\right)}{c}_1^3{\text{e}}^{3t}\\ y_{out}=\epsilon c_1{\text{e}}^t-\frac{a+1}{6\left(a+4\right)}{c}_1^3{\text{e}}^{3t}\\ z_{out}=\frac{1}{3}{c}_1^2{\text{e}}^{-t}{\text{e}}^{2t}-\frac{1}{15}\frac{a+1}{a+4}{c}_1^4{\text{e}}^{4t}\end{array}$ (22)

where $c_1$ can be user defined. Similarly, the expressions of the incoming vector can be set, as well. There, the integration constants associated with the “unstable eigenvalues” must be set equal to zero. The effectiveness of high order BC implemented compared to the classic first order ones, often encountered in bibliography, is presented in Figure 2 (approximation of outgoing solution) and in Figure 3(a) and Figure 3(b) (approximation of incoming solution).

Via the method of orthogonal collocation on finite elements combined with the aforementioned fourth order BC, the homoclinic connecting orbit of interest (i.e. at the origin) has been computed inside the truncated time interval $\left[-23.5,23.5\right]$ , which has been determined with the aid of the well-known Beyn’s method [16] , while the Euclidean distance of the initial point (i.e. the one closest to the equilibrium) of the outgoing solution vector from the equilibrium is ${\delta }_u\simeq 2.5545\times {10}^{-11}$ and the corresponding distance of the incoming (stable) solution vector is ${\delta }_s\simeq 4.2672\times {10}^{-10}$ , for the value of the active parameter $a\simeq \text{1}\text{.724323}\cdots$ . The trajectory of the homoclinic orbit is presented in Figure 4 together with the orbit resulting after numerical integration by use of the standard algorithm of numerical integration predictor-corrector method Adams-Bashforth-Moulton (ode113 of Mathworks Matlab).

3. Conclusion-Discussion

An efficient custom algorithm of orthogonal collocation on finite elements implemented in Mathworks Matlab has been applied to a laser model based on a modified version of Shimizu-Morioka system for the numerical location of a homoclinic orbit. An initial approximation of this orbit has become available through a numerical continuation of limit cycles, bifurcated from a Hopf bifurcation, up to a high value of the fundamental period. The efficiency of the algorithm lies in the fact that all the required equations, be them the collocation equations, the continuity equations, the boundary conditions and the phase condition, are converted to Matlab functions automatically, so that integrated, sophisticated Matlab routines used for solving systems of nonlinear algebraic equations, as well as optimization routines or any other relevant, user-defined routines can be applied directly for the solution of the aforementioned system of nonlinear algebraic equations. Furthermore, the high order boundary conditions defined and used herein can be useful when ordinary PCs of low to moderate computational power are used for the location of homoclinic orbits, as they do not require the increase of the length of the truncation interval in order to achieve the precision sought for the computation. The effectiveness of the high

order BC is evident in Figure 2, Figure 3(a) and Figure 3(b), where a backward (Figure 2) or forward (Figure 3(a) and Figure 3(b)) integration of the corresponding highest order approximation approaches the fixed point associated with the connecting orbit of interest satisfactorily.

Last, the physical meaning of the homoclinic orbit of the laser computed above has already been mentioned and it concerns the description of the saturation phenomenon present in this laser type. More precisely, the homoclinic point-to-point connecting orbit computed can be considered as a stylized representation of a high power self-pulsation associated with the aforementioned phenomenon (see Figure 4). Moreover, this application serves as a validation study of the implemented algorithm and the whole methodology on which the numerical procedure is based, as well.

Acknowledgements

P. S. Douris is pleased to acknowledge his financial support from “Andreas Mentzelopoulos Scholarships for the University of Patras”.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Cite this paper

Douris, P.S. and Markakis, M.P. (2018) Computation of a Point-to-Point Homoclinic Orbit for a Semiconductor Laser Model. Applied Mathematics, 9, 1258-1269. https://doi.org/10.4236/am.2018.911082

References