A New Chaotic System with Multiple Attractors: Dynamic Analysis, Circuit Realization and S-Box Design

This paper reports about a novel three-dimensional chaotic system with three nonlinearities. The system has one stable equilibrium, two stable equilibria and one saddle node, two saddle foci and one saddle node for different parameters. One salient feature of this novel system is its multiple attractors caused by different initial values. With the change of parameters, the system experiences mono-stability, bi-stability, mono-periodicity, bi-periodicity, one strange attractor, and two coexisting strange attractors. The complex dynamic behaviors of the system are revealed by analyzing the corresponding equilibria and using the numerical simulation method. In addition, an electronic circuit is given for implementing the chaotic attractors of the system. Using the new chaotic system, an S-Box is developed for cryptographic operations. Moreover, we test the performance of this produced S-Box and compare it to the existing S-Box studies.


Introduction
The discovery of the well-known Lorenz attractor [1] in 1963 opened the upsurge of chaos research. In the decades thereafter, a large number of meaningful achievements on chaos control, chaotification, synchronization and chaos application have emerged continuously. Great changes have also been made to the understanding of chaos. Scholars began to think more about a way to produce chaos rather than blindly suppress chaos. The generation of chaotic attractors in three-dimensional autonomous ordinary differential systems has been of particular interest. As we all know, a multitude of typical systems with chaotic attractors were found, including Rössler system, Chen system, Sprott system, Lü system, etc. [2][3][4][5][6][7][8].
With the further research of chaos, scientists found that some nonlinear dynamic systems not only have a chaotic attractor but also coexist with multiple attractors for a set of fixed parameter values. The coexisting attractors may be fixed points, limit cycles, strange attractors, etc. The number and type of attractors are usually associated with parameters and initial conditions of the system. Each attractor has its own basin of attraction which is composed of the initial conditions leading to long-term behavior that settles onto the attractor. The phenomenon of multiple attractors can be seen in many biological systems and physical systems [9][10][11]. In recent years, the low-dimensional autonomous chaotic systems with multiple attractors have aroused scholars' research enthusiasm.
with state vector (x, y, z) ∈ R 3 and parameter vector (a, b, c, k) ∈ R 4 . A butterfly attractor can be observed by numerical simulation on Matlab software (Matlab 8.0, MathWorks, Natick, MA, USA). The phase portraits of system (1) under parameters (a, b, c, k) = (4, 9,4,4) and initial condition (1, 1, 1) are shown in Figure 1. It visually demonstrates that system (1) displays an attractor as the system trajectories will eventually move to a bounded region. The Lyapunov exponents of the system are calculated as l 1 = 1.7729, l 2 = 0.0000, l 3 = −7.5949. The Lyapunov dimension is D l = 2 − l 1 /l 3 = 2.2334, so it can be determined that the attractor is a chaotic attractor. The time series of z generated from two very close initial conditions (1, 1, 1) and (1, 1, 1.001) are plotted in Figure 2. At the beginning, they are almost the same, but their differences are increasing after a number of iterations. That is to say, system (1) is sensitive dependence on initial conditions and its future behavior is unpredictable in the long term. The Poincaré map of system (1) is obtained via selecting the sections ∆ 1 = {(x, y) ∈ R 2 |z = 10} and ∆ 2 = {(y, z) ∈ R 2 |x = 0} . As shown in Figure 3, the Poincaré map is a sheet of point set. It is consistent with the nature of chaos.

The Stability of Equilibria
Suppose that parameters a, b, c, k are all positive real numbers. The equilibria of system (1) can be obtained by solvingẋ =ẏ =ż = 0. If k ≥ c √ ab, system (1) has only one equilibrium O(0, 0, k/c). If k < c √ ab, system (1) has three equilibria as follows:

Proposition 1.
Suppose that b > a > 0, k > 0, and the parameter c satisfies the following condition: then the equilibria O 1 and O 2 of system (1) are asymptotically stable.
Proof. By linearizing the system (1) at the equilibrium, the Jacobian matrix is given by By using |λI − H| = 0, the corresponding characteristic equation evaluated at O 1 , O 2 is obtained as where According to the Routh-Hurwitz criterion, the equilibria O 1 , O 2 are stable if all the roots of Equation (4) have negative real parts. This requires that w 1 > 0, w 2 > 0, w 3 > 0 and w 1 w 2 > w 3 . It is easy to verify that w 1 > 0, To make w 1 w 2 > w 3 , the parameter c should meet Since When the parameter c passes through the critical value c 0 , then double Hopf bifurcation occur with two limit cycles branched from O 1 , O 2 and system (1) loses its stability.
If c > k √ ab , then Equation (7) has a root with a positive real part.
all the roots of Equation (7) have negative real parts, which implies that O is stable.
It can be verified that O is asymptotically stable by applying the center manifold theorem.

The Evolution of Multiple Attractors
Detailed investigation of the complex dynamic behaviors of system (1) is presented in this section. Simulation experiments including bifurcation diagrams, phase portraits, Lyapunov exponents, and Poincaré maps give a close and intuitive look at system (1). There is a wealth of chaotic dynamics associated with the fractal properties of the attractor in system (1). With the change of parameters, system (1) experiences stable state, periodic state and chaotic state. For different initial values, system (1) performs different types of attractors with independent domains of attraction.

Dynamic Evolution with Parameter c
Consider the dynamic evolution of system (1) with respect to parameter c under the given parameter conditions a = 2, b = 8, k = 4. The bifurcation diagrams of system (1) versus c ∈ (0, 6) are shown in Figure 4a, where the red color branch and blue color branch are yielded from initial values x 01 = (1, 1, 1), x 02 = (−1, −1, 1), respectively. The overlapped regions of the red color and blue color branches indicate that the trajectories of x 01 , x 02 eventually tend to the same attractor, while the separated regions indicate that the trajectories of x 01 , x 02 tend to different attractors. Figure 4b is the Lyapunov exponents of system (1) with initial value x 01 . It shows that the system (1) experiences stable state, periodic state, chaotic state with the variation of c. When c ∈ (0, 1), system (1) is mono-stable as it has only one stable equilibrium. When c ∈ (1, 1.396), system (1) performs bi-stability with respect to the existence of two stable equilibria. As c increases across the critical value c 0 = 1.396, system (1) occurs double Hopf bifurcation at the equilibria. When c ∈ (1.396, 1.516), system (1) performs bi-periodicity. When c ∈ (1.516, 2.257), system (1) changes into mono-periodic state. When c ∈ (2.821, 3.096), system (1) yields two strange attractors from initial values x 01 , x 02 . When c ∈ (3.310, 4.167) ∪ (4.370, 5.324) ∪ (5.480, 6), system (1) has only one chaotic attractor. Table 1 describes the attractors of system (1) with different values of c. The phase portraits in Figure 5 illustrate the existence of different types of attractors in system (1).  (1)

Dynamic Evolution with Parameter k
The bifurcation diagrams of system (1) for parameters (a, b, c) = (4, 9, 4), k ∈ (5, 25) are shown in Figure 6a, where the red color branch and blue color branch are yielded from initial values x 01 , x 02 , respectively. Obviously, the state of system (1) changes from chaos to period and then to stable when parameter k increases from 5 to 25. It also can be illustrated by the Lyapunov exponents in Figure 6b

Electronic Circuit Realization
There are many works that are related to chaos based applications in the literature [31][32][33][34][35][36]. Here, we will present the circuit realization of system (1) for realistically obtaining its chaotic attractors. The numerical simulation in Figure 5e displays two coexisting strange attractors in system (1) for (a, b, c, k) = (2, 8, 2.9, 4) and initial conditions (±1, ±1, 1). For circuit realization of this state of system (1), we need to refrain from saturation of circuit elements, and the effective way to achieve this goal is to reduce the voltage values of the circuit via scaling the variables of system (1). In the process of scaling, we assume (X, Y, Z) = (x, y, z/2) and then the scaled system is obtained as (8) Figure 9 gives the new phase portraits of the scaled system (8) for (a, b, c, k) = (2, 8, 2.9, 4). Evidently, the scaling process does not cause fundamental changes to the system (1), but just limits the variables to a smaller region Ω = {(x, y, z) |x, y ∈ (−5, 5), z ∈ (0, 10)}.
The circuit diagram of system (8) raised by the OrCAD-PSpice programme (OrCAD 16.6, OrCAD company, Hillsboro, OR, USA) is presented in Figure 10. It has three input (or output) signals with respect to the variables X, Y, Z, and the operations between signals realized via the basic electronic materials including resistors, capacitor, TL081 operational amplifiers (op-amps), and AD633 multipliers. By fixing R1 = R3 = 20 KΩ, R2 = 200 KΩ, R4 = 50 KΩ, R5 = R6 = 100 KΩ, R7 = 138 KΩ, R8 = 3000 KΩ, R9 = 4 KΩ, C1 = C2 = C3 = 1 nF, Vn = −15 V, V p = 15 V and executing the circuit on electronic card shown in Figure 11, we can obtain the outputs of circuit in the oscilloscope. The oscilloscope graphics in Figure 12 show good consistency with the numerical simulations in Figure 9. Hence, we can come to a conclusion that the coexisting chaotic attractors in system (1) are physically obtained.

S-Box Design and Its Performance Analysis
This section aims to raise a new chaotic S-Box algorithm by applying system (1). In the algorithm design, first random number generation is performed and then S-Box is produced. The S-Box generation algorithm pseudo code is shown in Algorithm 1. For establishing the S-Box production algorithm, we first input the parameters (a, b, c, k) = (2, 8, 2.9, 4) and initial value (x 0 , y 0 , z 0 ) = (−1, −1, 1) of the system and then the float number outputs are produced. In order to generate more random outputs in the analysis of the chaotic system, we select an appropriate step interval ∆h and used it as a sample value. More random sequences are obtained by setting the appropriate step interval ∆h = 0.000001. System (1) is solved by using the RK4 algorithm with the initial conditions and the specified sampling value, and time series are obtained. In our designed chaotic S-Box algorithm, the outputs of y, z phases of system (1) are used. Float number values (32 bits) obtained from these phases are converted to a binary system. By taking 8 bits from the low significance parts (LSB) of the 32-bit number sequences generated from both phases, these values are XORed. The obtained new 8-bit value is converted to a decimal number. This value is discarded if the decimal number was previously generated and included in the S-Box, if not produced before, it is added to the S-Box. In this way, this process continues until the distinct 256 values (between 0 and 255) are obtained on the S-Box. The generated S-Box is shown in Table 2.

Algorithm 1
The S-Box generation algorithm pseudo code. 1: Start; 2: Inputting parameters and initial value of the system; 3: Sampling with step interval ∆h; 4: i = 1, S-Box = []; 5: while (i < 257) do 6: Solving system with RK4 algorithm and obtaining time series (y, z); 7: Convert float to binary number; 8: Take LSB-8 bit value from RNG y ⊕ z phase; 9: Convert binary to decimal number (8 bit) 10: if (Is there decimal value in S-Box = yes) then 11: Go step 6.  In order to determine that the produced S-Boxes are robust and strong against attack, some performance tests are applied. We mainly focus on these tests: nonlinearity, outputs' bit independence criterion (BIC), strict avalanche criterion (SAC), and differential approach probability (DP). In addition, the comparisons of the performance between this new S-Box and the existing chaotic S-Box proposed by Chen [37], Khan [38], Wang [39], Ozkaynak [40], Jakimoski [41], Hussain [42], Tang [43] are presented in Table 3.
Nonlinearity is regarded to be the most core part of all the performance tests. The nonlinearity values of the S-Box yielded by system (1) are obtained as 104, 106, 104, 104, 108, 104, 110 and 104. Accordingly, its average value, minimum value and maximum value are computed as 105, 104 and 110.
By comparing the nonlinearity of other S-Box shown in Table 3, we can claim that the new S-Box is better than others in some measure. Table 3. The comparison of different chaotic S-Boxes (BIC: bit independence criterion; SAC: strict avalanche criterion; DP: differential approach probability).  [44]. Generally speaking, the establishment of SAC implies a possibility that half of each output bit will be changed with the change of a single bit. Table 3 tells the average, minimum and maximum SAC values of the new S-Box as 0.5014, 0.3906, 0.5937. Evidently, the average value of the new S-Box is close to the ideal value 0.5. BIC is also an important criterion found by Webster et al. [44]. It can partially measure the security of cryptosystems. The set of vectors generated by reversing one bit of the open text is tested to be independent of all the pairs of avalanche variables. While the relation between avalanches is measured, variable pairs are necessary to calculate the correlation value [45]. BIC-SAC and BIC-Nonlinearity values are calculated when the BIC value is calculated. When the values in Table 3 are examined, the BIC-SAC values are calculated as follows: average value 0.5028, minimum value 0.4394 and maximum value 0.5312. The average value is almost equal to the optimum value 0.5. DP is another performance index for testing the S-Box, which is established by Biham et al. [46]. In this analysis, the XOR distribution balance between the input and output bits of the S-Box is determined. The very close probability of XOR distribution between input and output bits often indicates the ability to resist the differential attack of the S-Box. The low DP value suggests that the S-Box is more resistant to attack. The minimum and maximum DP values of the new S-Box are determined as 4.0 and 10. From Table 3, we know that the DP value of the new S-Box is the same as the S-Boxes presented by Tang and Ozkaynak. After testing the performance of the new S-Box by using some important indices and comparing with other S-Boxes, we can determine that the new S-Box generated by system (1) has better performance than other S-Boxes. Thus, it will be more suitable for attack resistant and strong encryption.

Conclusions
A special chaotic system with multiple attractors was studied in this letter. The complex dynamic behaviors of the system were mainly presented by numerical simulations. Bifurcation diagrams and phase portraits indicated that the system exhibits a pair of point attractors, a pair of periodic attractors, and a pair of strange attractors with the variation of system parameters. In addition, an electronic circuit was designed for realizing the chaotic attractors of the system. Moreover, a new S-Box was generated by applying the chaotic system, and the performance evaluation and comparison of the S-Box were presented. It showed that the new S-Box has better performance than some existing S-Boxes. Actually, the study of chaotic system with multiple attractors is of recent interest. More important issues corresponding to this topic will be addressed in our future paper.