Dynamics and Complexity of a New 4D Chaotic Laser System

Derived from Lorenz-Haken equations, this paper presents a new 4D chaotic laser system with three equilibria and only two quadratic nonlinearities. Dynamics analysis, including stability of symmetric equilibria and the existence of coexisting multiple Hopf bifurcations on these equilibria, are investigated, and the complex coexisting behaviors of two and three attractors of stable point and chaotic are numerically revealed. Moreover, a conducted research on the complexity of the laser system reveals that the complexity of the system time series can locate and determine the parameters and initial values that show coexisting attractors. To investigate how much a chaotic system with multistability behavior is suitable for cryptographic applications, we generate a pseudo-random number generator (PRNG) based on the complexity results of the laser system. The randomness test results show that the generated PRNG from the multistability regions fail to pass most of the statistical tests.


Inroduction
The chaotic behavior as a rich nonlinear phenomenon has been detected in many non-natural and natural systems, and usually plays an important role in their performance [1,2]. Chaotic systems are complicated and have many interesting features, such as unpredictability, topological mixing, and high sensitivity to their initial conditions and parameters [3,4]. Therefore, chaotic systems have received significant attention from various fields including cryptography [5,6], secure communications [7,8], laser applications [9,10], biomedical engineering [11,12], and many others.
Existing chaotic systems can be classified into two categories: systems with self-excited attractors and systems with hidden attractors [13]. The chaotic system with self-excited attractors has a basin of attraction that is intersected with an unstable equilibrium, whereas the chaotic system with hidden attractors has a basin of attraction which does not intersect with any open neighborhoods of equilibria [14,15]. According to the above definition, most of the classical chaotic attractors are self-excited [16,17]. Meanwhile, it has been demonstrated that the attractors in dynamical systems with no equilibria [18,19], stable equilibria [20], lines of equilibria [21], and curves of equilibria [22] are hidden attractors.
However, with further investigation of chaos, it was unexpected to find that many systems with self-excited and hidden attractors have more than one attractor for a given set of parameters and (i) We derive a new 4D chaotic laser system with three equilibria from Lorenz-Haken equations; (ii) We investigate the stability of the symmetric equilibria, and the existence of coexisting multiple Hopf bifurcations on these equilibria; (iii) We analyze the presence of complex coexisting behaviors in the laser system; (iv) We use the complexity of the laser system time series to locate the regions of coexisting attractors when the parameters and initial values vary; (v) Based on the complexity of the system time series, we study the randomness of multistability regions.
The rest of this paper is organized as follows: Section 2 introduces the new 4D chaotic laser system and studies its dynamical properties. Section 3 investigates the existence of Hopf bifurcation in the laser system. Section 4 provides the details about the multistability of the laser system. In Section 5, we use SamEn to locate the regions of the coexisting attractors, as well as to demonstrate the randomness of these regions. The conclusions are presented in Section 6.

A New 4D Chaotic Laser System From Lorenz-Haken Model
In this section, we discuss the dynamics of a new 4D chaotic laser system which is derived from the well-known Lorenz-Haken equations [35]. In the standard notation of Reference [36], the Lorenz-Haken equations is given by In the optical language, x is proportional to the electric field, y is proportional to the induced macroscopic polarization, (r − z) denotes the inversion, σ = τ P τ E , and b = τ P τ N . Here, τ E represents the optical field, τ P is the induced polarization, and τ N denotes the inversion parameter. Meanwhile, the parameter δ governs the coupling between amplitude and phase variations, and q is known as the linewidth enhancement factor.
Since both x and z can be chosen to be real [37], the dynamics of Equation (1) can be investigated by considering the following linear transformation Consequently, the new 4D chaotic laser system is defined as where x i are state variables and σ, δ, r and b are parameters.

Dissipation and Symmetry
The divergence of system (2) is defined as Thus, the system (2) becomes dissipative when (σ + b + 2) > 0. This means each volume element V 0 e −(σ+b+2)t of system (2) shrinks to zero as t −→ ∞ at an exponential rate (σ + b + 2). Additionally, the system (2) has invariance under the coordinate transformation Consequently, the system (2) has rotational symmetry around the x 4 -axis.

Equilibria and Stability
Suppose that the parameters σ > 0, δ > 0, r > 0 and b > 0, then the equilibria of the system (2) can be calculated by solving the following equations From the above equations, it can be obtained that the equilibria of the system (2) have the following form: where k is either 0 or ± b(r − (1 + δ 2 )). The system (2) has one real equilibrium E 1 (0, 0, 0, 0) when r = 1 + δ 2 , whereas it has three real equilibria if r > 1 + δ 2 Using the Jacobian matrix, the system (2) is linearized at the equilibrium E i as follows Since the equilibria E 2,3 are symmetric about the x 4 −axis, then they will have the same characteristics. Therefore, the characteristic equation of Jacobian matrix at the equilibrium E 2 with b = 1 can be written as where It is obvious that Equation (3) always has one eigenvalue with negative real part which is λ 1 = −1, whereas the real parts of the other eigenvalues are not always negative. It is well-known that a system is asymptotically stable when all eigenvalues have negative real parts; otherwise, the system is unstable. By Routh-Hurwitz criterion, the real parts of all the eigenvalues of the system (2) are negative if and only if By choosing the parameters σ > 0 and δ > 0, these inequalities lead to the following condition: Thus, if the above conditions are satisfied, then the equilibrium E 2 is an asymptotically stable.

Local Bifurcation Analysis and Numerical Simulations
This section reviews the Hopf bifurcation using the bifurcation theories. In addition, the existence of coexisting symmetric Hopf bifurcations in the system (2) will be investigated with the variation of parameter r ∈ R + .

Hopf Bifurcation
Hopf bifurcation is the source of a limit cycle, which usually appears when the stability of the equilibrium point changes at some critical parameter value. To illustrate the Hopf bifurcation of a dynamical system on the equilibrium point, consider a vector field as followṡ where x ∈ R 4 and ζ ∈ R + represent the phase variables and the parameters, respectively. The vector field undergoes a Hopf bifurcation when the following conditions are satisfied simultaneously [38]: (A) nondegeneracy condition: the Jacobian matrix J (x 0 ,ζ 0 ) has one pair of purely imaginary roots, and other roots have nonzero real parts; (B) transversality condition: the real part of differentiation characteristic equation with respect to the parameter ζ satisfy (C) the first Lyapunov coefficient l 1 is nonzero.
In order to derive the first Lyapunov coefficient l 1 , suppose that Equation (2) has an equilibrium point at x = x 0 . By denoting X = x − x 0 , we can write as where A is the Jacobian matrix, and B and C are symmetric multilinear vector functions which are defined as Suppose that A possesses a pair of purely imaginary eigenvalues λ 1,2 = ±iω, meanwhile, the other eigenvalues have nonzero real part. Let p, q be an eigenvectors of A satisfying the following three conditions By means of an immersion of the form X = V(µ, µ), the 2D center manifold associated to the eigenvalues λ 1,2 = ±iω is parameterized, where V : C 2 −→ R 4 has a Taylor expansion of the following form with v jk ∈ C 4 and v jk = v jk . By substituting Equation (12) into (8) Defined by the coefficients µ j µ k , the complex vectors v jk can be obtained by solving Equation (13). On the chart µ for a center manifold, the system (13) can be written aṡ Thus, the first Lyapunov coefficient can be defined as where

Numerical Simulations
To investigate the existence of Hopf bifurcation in the system (2) at the equilibrium E 2 , we will examine the conditions (A), (B) and (C) one by one.
At last, we will calculate the first Lyapunov coefficient l 1 under the above fixed parameters. The Jacobian matrix J on the equilibrium point E 2 is given by The proper eigenvectors q and p are obtained by straightforward calculations where the above eigenvectors q and p satisfy the three conditions (11), namely From Equation (10), the multilinear vector functions of the system (2) are calculated as follows From (22)- (24), it follows that Thus, one obtains Consequently, the first Lyapunov is obtained by substituting (26) into (15) Therefore, the Hopf bifurcation of the system (2) at equilibrium point E 2 is nondegenerate and supercritical. Furthermore, the equilibria E 2 and E 3 are symmetric about the x 4 −axis, hence, the system (2) should also undergo a Hopf bifurcation at E 3 . Two numerical simulations are given in Figure 3. For r = 5.5 < r 0 , the orbit of the system (2) with the initial values (1.8, 1.8, 2, 4) is attracted to the stable equilibrium point E 2 , whereas the orbit with the initial values (−1.8, −1.8, −2, 4) is attracted to the other stable equilibrium point E 3 , as illustrated in Figure 3a. In Figure 3b, by choosing r = 6.5 > r 0 with the initial values (1.8, 1.8, 2, 4) and (−1.8, −1.8, −2, 4), the orbits of the system are attracted to stable limit cycles emerging from E 2 and E 3 , respectively.
According to Reference [39], m = 2, τ = 1 and r = 0.1 ∼ 0.2 times standard deviation (SD) of the time series. In our experiment, we fix m = 2, τ = 1 and r = 0.2 × SD.  (2): (a) r = 5.5 < r 0 , the orbit of the system is attracted to the stable symmetric equilibria E 2 and E 3 ; (b) r = 6.5 > r 0 , the orbit of the system is attracted to a stable limit cycle emerging from the symmetric equilibria E 2 and E 3 .

Multistability Behavior
A nonlinear dynamical system with multistability behavior can generate two or more attractors simultaneously depending on the initial values of the system. This section investigates the existence of multistability behavior in the system (2).
When we fix the parameters σ = 2, δ = 1.5, b = 0.7 and select r as bifurcation parameter for over the range r ∈ [7.5, 10], the coexisting bifurcation models of the state variable x 1 are depicted in Figure 4a. In this figure, the attractor colored in blue is initiated from (−2, 1, 1, 1), meanwhile the attractor colored in red begins with the initial conditions (1, 1, 1, 1). As can be observed in Figure 4a, the system (2) shows coexisting multiple chaotic attractors as well as the coexistence of multiple quasi-periodic attractors. To show the coexistence of multiple chaotic attractors visually, Figure 5 plots different orientations of the phase portraits of the system (2) when its parameters are set as σ = 2, δ = 1.5, b = 0.7, and r = 9.41.  In addition, when we set σ = 4, δ = 0.5, b = 2 with 26 ≤ r ≤ 30, Figure 4b shows that the chaotic attractor with two stable fixed-point attractors coexist for the initial values (±2, 1, 1, ±2). For the orbit colored in blue, the evolution begins from attracting to the stable equilibrium E 3 within the range 26 ≤ r ≤ 26.7, and then the system shows chaotic behavior when r ≥ 26.8. For (−2, 1, 1, −2) (red), the system converges to the stable equilibrium E 2 when 26 ≤ r ≤ 28, and then exhibits chaotic behavior when r ≥ 28.1. For the initial values (2, 1, 1, −2) (green), the system attracts to the stable equilibrium E 3 when 26 ≤ r ≤ 27.8, meanwhile the chaotic behavior is shown when r ≥ 27.9. Selecting r = 27, an interesting dynamic is observed in the system (2) by plotting different orientations of the phase portraits with the corresponding time series, as shown in Figure 6. These portraits confirm the coexistence of three different attractors: (a) blue butterfly attractors surrounds the symmetric equilibria E 2 and E 3 ; (b) the red stable fixed-point attractor for E 2 , and the green stable fixed-point attractor for E 3 .
Through the above analysis, we can observe that the multistability behavior occurs in the system (2) with various kinds of coexisting attractors. Therefore, it can be concluded that the system (2) has high sensitivity to both initial values and parameters.

Complexity and Randomness of Multistability Regions
This section discusses determining and locating the parameters and initial values that show multistability behaviors, as well as investigates the randomness of the multistability regions.

Sample Entropy
Sample Entropy (SamEn) is a mathematical algorithm proposed by Richman [40]. It is used to provide a quantitative explanation about the complexity of nonlinear dynamical systems. Obviously, a system with bigger SamEn values indicates that it requires additional information to predict its attractor, hence, it is a chaotic system. Suppose that the time series (y i , i = 0, 1, 2, . . . , M − 1) of a dynamical system with a length of M, then the SamEn algorithm can be calculated by the following steps: (A) Reconstructing phase-space: for a given embedding dimension m and time delay τ, the reconstruction sequences are given by where i = 1, 2, . . . , M − m + τ. (B) Counting the vector pairs: let B i be the number of vector Y j such that where r is the tolerance parameter, and d[Y i , Y j ] is the distance between Y i and Y j , which is defined by (C) Calculating probability: according to the obtained number of vector pairs, we can obtain then calculate the probability by (D) Calculating SamEn: repeating the above steps we can obtain φ m+1 (r), then SamEn is given by According to Reference [39], m = 2, τ = 1 and r = 0.1 ∼ 0.2 times standard deviation (SD) of the time series. In our experiment, we fix m = 2, τ = 1 and r = 0.2 × SD.
It is well-known that the cross-section of the basins of attraction can determine the dynamical system behaviors when its initial values vary. However, it is interesting to ask if there is any technique that can determine the behaviors of a dynamical system when its initial values and parameters vary. Therefore, SamEn based contour plots are applied to locate the regions of chaotic and periodic state, and hence, to determine the parameters and initial values that show multistability behaviors. To locate those parameters and initial values in the system (2), we designed the following experiments: (1) consider r as bifurcation parameter and set σ = 4, b = 2 and δ = 0.5; (2) let (x 10 , x 20 , x 30 , x 40 ) be the initial values; (3) calculate SamEn versus varying the parameter r and one of an initial value; (4) calculate SamEn versus varying two of the initial values. Figure 7 plots SamEn of the system (2) in a two-dimensional plane when r ∈ (24, 30) and different initial values. It can be observed from Figure 7a-d that four cases are analyzed when the initial values are set as (x 10 , 1, 1, 2), (2, x 20 , 1, 2), (2, 1, x 30 , 2) and (2, 1, 1, x 40 ), respectively. From Figure 7, it can be seen that the parameter r and the initial values in the blue regions have smaller SamEn values, which means that the system (2) shows periodic state, whereas, those in the yellow and green regions lead to a chaotic state. Furthermore, Figure 8 shows the chaotic and periodic regions of system (2) when two of the initial values vary simultaneously.

Chaos-Based PRNG
Many chaotic systems have been applied to generate pseudorandom number generator (PRNG). The need of PRNG arises in many cryptographic applications, e.g., common cryptosystems employ keys, data hiding, and auxiliary quantities used in generating digital signatures [41,42]. However, secret keys of most chaos-based cryptographic schemes are generated by parameters and initial values of the employed chaotic systems [43]. Those parameters and initial values might be from multistability regions; it is therefore important to investigate the randomness of the trajectories generating from multistability regions.
To investigate the randomness of blue-green regions (multistability behaviors) and green regions (chaotic), which is shown in Figure 7d, we use here a simple chaos-based PRNG as an example.
The generation procedures of the chaos-based PRNG are shown in Algorithm 1, for which x 1 , x 2 , x 3 and x 4 generates 1,000,000 bits binary string.
Several statistical tests can be employed to test the randomness of PRNG. Our experiment uses the highest standards of statistical packages which is NIST-800-22 [42]. The NIST-800-22 consists of 16 empirical statistical tests that provide true evaluation for the randomness of PRNG. Each test is developed to detect the non-random areas of a binary sequence from different sides, and then to derive a p-value. According to the recommendations in [24,44], we set the confidence level α = 0.01, and we use a binary sequence with length of 1,000,000 bit as the testing input. Since the confidence level of each test in NIST is set to be 1%, then the sequence is considered to be random with a confidence of 99% when the obtained p-value is bigger than 0.01.
According to Algorithm 1, we can obtain four PRNG from the trajectory of x 1 , x 2 , x 3 and x 4 when the initial values are considered as input. For σ = 4, δ = 0.5, b = 2 and r ∈ [27,29] with the initial values (2, 1, 1, −2), the SamEn values of the selected parameters and initial values are within the blue-green regions (multistability), as shown in Figure 7d. The randomness of the corresponding PRNG that generated from the trajectory of x 1 , x 2 , x 3 and x 4 can be visually shown by depicting the NIST-800-22 test results, as seen in Figure 9. As can be observed from Figure 9, the four PRNG generating from multistability regions fail to pass most of the statistical tests. On the other hand, when σ = 4, δ = 0.5, b = 2 and r ∈ [27,29] with the initial values (2, 1, 1, 2), the SamEn values are within the green region (chaotic), as shown in Figure 7d. Table 1 lists the corresponding NIST-800-22 results for each of the four PRNG. It is obvious that the four PRNG can pass all the statistical tests. ApEn cumul. f. cumul. r. lempel-ziv non-overlap universal runs binary rank overlap Figure 9. The statistical tests NIST SP800-22 of the pseudorandom number generator (PRNG) that generated by x 1 , x 2 , x 3 , x 4 of the system (2) with σ = 4, δ = 0.5, b = 2, r ∈ [27,29] and for the initial values (2, 1, 1, −2).

Conclusions
This paper has introduced a new 4D chaotic laser system, which is derived from Lorenz-Haken equations. The new chaotic laser system has three equilibria and only two quadratic nonlinearities. The dynamics of the new system have been studied deeply, in which the system shows coexisting multiple Hopf bifurcations, and complex coexisting behaviors of two and three attractors. In addition, we applied SamEn contour plots for measuring the complexity of the system when its initial values and parameters vary. Simulation results have shown that multistability regions can be easily determined and located using SamEn contour plots. To examine the randomness of PRNG that generate from the multistability regions, we used the NIST-800-22 tests. Statistical test results demonstrate that the generated PRNG from multistability regions are non-random. This means that although the multistability behaviors indicate high sensitivity of chaotic systems, they might be unsuitable for cryptographic applications.