Generating Infinitely Many Coexisting Attractors via a New 3D Cosine System and Its Application in Image Encryption

In this paper, a new third-order chaotic system which has extremely multistability is constructed by introducing the boosted control of cosine function. In comparison with other chaotic systems of extremely multistability, the proposed chaotic system can spontaneously generate the infinitely many coexisting attractors towards two directions of the phase plane. It indicates the proposed system can output more chaotic sequences of different amplitudes at the same time. This peculiar physical phenomena is very interesting and worth studying. Relative to original chaotic system, the chaos characteristic of the proposed system is obviously enhanced, the value of max Lyapunov exponent is increased significantly and the complexity value was higher. In particular, many periodic windows of the original chaotic system become chaos. It means the proposed chaotic system has better chaotic characteristics. If the new system is applied to the field of cryptography, it would be a better system model as a pseudo-random signal generator (PRSG). Then, the new image encryption algorithm is designed based on the proposed discrete system, and its safety performance is tested. The experimental results demonstrate the feasibility of its application in the field of cryptography.


I. INTRODUCTION
Chaoic system has received a lot of attentions on account of its bright application prospects in the field of nonlinear engineering [1]- [6]. Generally speaking, extremely multistability shows that there are infinitely many numerical solution in the differential equations of a dynamical system with varying initial states [7], [8]. The original chaotic system of infinitely many equilibrium points or switchable equilibrium point has the ability to generate infinitely many coexisting attractors [9], [10]. Recently, another better method for constructing the system with extremely multistability was found, the initial-offset boosted attractors of infinitely many coexistence can be obtained by introducing the trigonometric functions to some specific linear terms of chaotic system [11]- [13]. The chaotic system which has infinitely many coexisting attractors depends extremely on the initial state. Thus, in comparison with some single stable The associate editor coordinating the review of this manuscript and approving it for publication was Bing Li . systems [14], [15], its dynamic behaviors are often more complex.
The initial-offset boosted multistability can provide more ergodicity and flexibility for some engineering applications based on chaos theory [16]. In particular, it can be applied to switch different dynamic states in the control of multistable state if combined with several suitable algorithms of control science [17]. Some recent studies show that a part of classic chaotic systems have the ability to generate coexisting attractors by their standalone attractor basins [18]- [20]. Classical 3D Lorenz system can produce two symmetric attractors from single attractor by adjusting the initial conditions [21]. Ye et al. found a new circuit system of extremely multistability by introducing a meminductive model to the 4D Wienbridge oscillator [22]. In particular, some chaotic systems can generate extremely multistability by the boosted control of trigonometric function. Lai et al. obtained infinitely many coexisting attractors by introducing sine function to the Sprott B system [23]. The discovery of extremely multistable system was a leap from the finite to the infinite. Its potential industrial application value is worth further and deeply exploring.
Third-order chaotic system is the lowest order continuous system which can generate the chaos attractor [24], it is the origin of research on continuous chaotic system. However, due to the lack of multiple feedback of the nonlinear term, 3D chaotic system is hard to spontaneously generate extremely multistability. Thus, constructing the extremely multistability based on 3D chaotic system is very crucial for chaotic basic demonstration teaching. This paper successfully constructed a new 3D chaotic system which can produce infinitely many coexisting attractors by using the boosted control of cosine function. Especially, due to the introduction of two cosine functions, it can generate two directions infinitely many coexisting attractors by changing different initial conditions. That means the proposed system can simultaneously output the pseudo-random signals of two directions in the phase plane. It greatly increased the chaos range and enhanced the pseudo-random performance of the chaos sequence. The proposed chaotic system has better chaos characteristic and higher sequence complexity, it provided a good model for the field of cryptography.
In this article, we focus on a new chaotic system of extremely multistability which can generate infinitely many coexisting attractors. It is organized as follows. In section 2, a new system model is proposed based on sprott system, and the equilibrium points is calculated. In section 3, Lyapunov exponents and Kaplan-Yorke dimension of the new system is calculated, its 3D chaotic attractor is presented. In section 4, the dynamic behavior of enhanced chaos in this proposed system is analyzed, and the complexity is compared. In section 5, the dynamic behavior of infinitely many coexisting attractors is analyzed. In section 6, the new encryption algorithm is designed, and its safety performance is tested. Finally, the proposed new system is successfully implemented by DSP.

II. CHAOTIC SYSTEM MODEL A. SYSTEM EQUATIONS
The proposed system is a third-dimensional chaotic system of extremely multistability, and it originates from the restruction for a 3D chaos system proposed by Ref. [24]. The mathematical model of Ref. [24] is: where, x, y, z are state variables, and a, b, c are system parameter. This 3D system cannot generate coexisting attractors by itself standalone attractor basins. Then, reconstructing the system equation by introducing two boosted control of cosine function for state variables x, z, the chaotic system can generate infinitely many coexisting attractors of two directions on the phase plane. The mathematical equations of the proposed system be expressed as: in which, x, y, z are all state variables, and a, b are the positive constants. Parameter a can control the amplitude of chaotic sequence and b as a parameter of the nonlinear feedback loop. In order to reduce the period of cosine function, expand state variables x, z as 2x, 2z. Its dynamic characteristics can be analyzed based on this chaotic system equations.

B. ANALYSIS OF DISSIPATION
If system (2) is a dissipation system, its condition is: Thus, when x ∈ (kπ + π/2, kπ + π) (k = 0, ±1, ±2. . . .), the new system is considered to be dissipative. It indicates the running track of the system will compress to an empty set, and the progressive motion will tends to be stable in a domain of attraction.

C. EQUILIBRIUM POINTS SET
Letẋ =ẏ =ż = 0, the equilibrium point of system (2) can be calculated: (4) in which, k = 0, ±1, ±2 . . . .., the Jacobian matrix at the equilibrium point S is: The third-order characteristic polynomial can be obtained: Thus, we can calculate: According to 3D Routh-Hurwitz stability criterion, we have the sufficient and necessary conditions for stability of the proposed new chaotic system is: Because the equilibrium points of the proposed chaotic system are changed with the periodic change of the cosine function. An infinite number of attractors may be generated due to the extremely multistability of the chaos system. VOLUME 9, 2021

III. ANALYSIS OF DYNAMIC BEHAVIOR A. LYAPUNOV EXPONENTS AND KAPLAN-YORKE DIMENSION
The sensitivity to initial conditions is an important dynamical characteristic of chaos. It can cause a larger separation of two near orbits with varying initial state. These characteristics can be described by Lyapunov exponents, the positive Lyapunov exponent shows the instability of the phase orbitals, it indicates the system is chaos. The negative Lyapunov exponent shows the contraction of the running orbitals, it reflects the unstable points or periodic points. The attractors are formed from the repeated folding of chaotic orbits. The mathematical formula of a continuous system of differential equations be mathematically defined as: where, Y (t) is a function with varying time,Ẏ (t) is a derivative of Y (t), F[Y (t)] means some functional relationship between two functions. Then, Eq. (9) can be treated by using Euler discretization method, and we have: t is a minute time variable. Define A p are the p-order compound matrices of Jacobian matrix generated by Y (t), then we have: in the last formula, A p (t) is a p-order compound matrices on time variable. K p [Y (t)] is the limit value of the p-order compound matrices of Jacobian matrix within the range of time variable will approach to 0. Finally, we can get the computational formula of Lyapunov exponent of the continuous system: where, S p (t) is the trace ofÃ p (t). The Lyapunov exponent value of the continuous chaos system can be calculated and obtained according to the above formula. Generally speaking, for a chaos system of n dimensions, the number of Lyapunov exponents are also n, the Lyapunov exponents λ i (i = 1,2. . . .n) need satisfy the conditions λ 1 > λ 2 > . . . . > λ n . If it satisfies: where, j is the max value of i which can satisfy the conditions λ i > 0, then we can obtain the Kaplan-Yorke dimension λ L of Lyapunov exponent by the following formula: Setting system parameter a = 15, b = 0.1, the initial conditions of system. (2) (x 0 , y 0 , z 0 )=(0,1,0), The Lyapunov Setting system parameter a = 15, b = 0.1, the initial conditions of system. (2) (x 0 , y 0 , z 0 )=(0,1,0). By using ODE45 algorithm, the chaotic attractors can be obtained by using numerical simulation. Fig. 1(a) shows the chaotic attractor in the 3D phase spaces. Obviously, the system is chaos in this case. To further explain the formation of attractors, the time sequences are analyzed and simulated. Let the time step is 0.01s and select the length 300s∼600s as the object of study. The result of numerical simulation is represented in Fig. 1(b), it shows the system is in a stable chaotic orbits, and the attractors are formed by the superposition of their chaotic oscillations.

IV. ENHANCED CHAOS OF BOOSTED CONTROL A. ENHANCED CHAOS IN THE DYNAMIC BEHAVIOR
Setting the initial value of system (1) (x 0 , y 0 , z 0 )=(0,1,0) and system parameter c = 1, b = 0.01. Parameter a is the variable parameter, and the varying range is [11,15]. As shown in Fig. 2(a), the bifurcation diagram of system (1) has a biggest periodic windows within a ∈ [12.735, 13.226]. Then, more periodic windows are generated as parameter a changes. When a ∈ [13.833, 14.198], the system has another larger window. The generation of these periodic windows causes the system to produce intermittent chaotic oscillations. It has a negative effect on the industrial application of chaotic system. Let the initial state of the proposed system (x 0 , y 0 , z 0 )=(0,1,0) and system parameter b = 0.01. Parameter a as the variable parameter, and the varying range is also [11,15].
(when system parameter of system (1) a = 1 and c = 0.01, the chaotic system of system (1) is equivalent to the proposed system without introducing two boosted control of cosine function) Fig. 2(b) shows the bifurcation diagram of the proposed chaotic system within the equivalent varying range of system (1). From the simulation result, most of periodic windows in Fig. 2(a) have become chaos due to the introduction of boosted control of cosine function. In comparison with system (1), the proposed chaotic system has almost no periodic windows within the test interval, and the chaotic state fills the whole range. That is to say, the chaotic behavior is enhanced. It greatly enhances the application of chaotic pseudo-random signals.
In order to further study the effect of system parameters on chaotic behavior. Keep the above parameters unchanged, select two system parameters a and b as the variable parameter, the varying range a ∈ [11,15] and b ∈ [0, 1] respectively. With a and b synchronously changed, the max Lyapunov exponent of system (1) is shown in Fig. 3(a). When a ∈ [11,15] and b ∈ [0.2, 0.4], system (1) has a bigger periodic windows. In particular, when b ∈ [0.6, 1], the max Lyapunov exponents of system (1) all go down to 0, it indicates the system is in a periodic state, and the max value of the max Lyapunov exponent in this test range is 0.332. Fig. 3(b) shows the max Lyapunov exponent of the proposed chaos system. For comparative analysis, select the same varying range a ∈ [11,15] and b ∈ [0, 1] as the test range. In comparison with system (1), the bigger periodic windows of a ∈ [11,15] and b ∈ [0.2, 0.4] has disappeared, the max Lyapunov exponent value in this test range are within 0.5∼1. It shows the periodic windows in this range has become chaos state. The dynamic behaviors of system (1) is periodic in the range b ∈ [0.6, 1]. However, by introducing two boosted control of cosine function to state variables x, z, the max Lyapunov exponent of the proposed chaotic system has been improved significantly. In this test interval, the max value of the max Lyapunov exponent of the proposed chaotic system is 1.198. Relative to system (1), the chaos interval of the proposed chaos system has also been significantly increased, the chaotic property is enhanced.

B. ENHANCED CHAOS IN THE COMPLEXITY
The enhancement of dynamic behavior makes chaotic behavior transition from periodic state to chaotic state. In terms of the application of chaotic system, the complexity characteristic of chaotic sequences can more directly reflects the usability of a chaos system. Spectral Entropy (SE) complexity algorithm is a measure of a system based on Shannon entropy, it can effectively shows the complexity characteristic of chaotic system. The higher value of complexity reflects the stronger chaotic behavior generated by the chaos system. Spectral Entropy (SE) complexity algorithm can be described:x where, x(n) is the energy amplitude of system. The function of Eq. (15) is to remove DC part of sequence. Then, do discrete Fourier transform, we have: in which, k = 0, 1, 2 . . . , N − 1. Then, by using the Parseval theorem, we can calculate the relative power spectrum: In Eq. (17), k = 0, 1, 2 . . . , N /2 − 1, and we have the total power:p Thus, the probability of power spectrum can be calculate, and we have: in which, N is the sequence length, X (k) is the sequence after FFT. Combined with the above formula, we can obtain the computational formula of SE: Based on the above formula, the spectral Entropy (SE) complexity of chaos sequence can be calculate, and the corresponding dynamic behavior can be analyzed.
C 0 complexity algorithm is also a complexity algorithm based on FFT. The proportion of irregular part between the regular and the irregular parts is C 0 complexity value. Its specific calculation method can be defined as: where, k = 0, 1, . . . , N − 1,X (k) is sequence the after FFT. Then, remove the regular part of sequence, we have: in which, Y N is mean square value. Then, if we introduce parameter r as a tolerance parameter, and reserve the mean square value which more than r, and other parts are zero, we have:X Then, do the inverse FFT, we can get: where, n = 0, 1, . . . , N − 1. According to the above equations, C 0 complexity can be obtain: In order to compare and analyze the above simulation results of bifurcation diagram and max Lyapunov exponent, keeping other parameters unchanged, select parameter a and b as the varying parameters, and the test ranges are also a ∈ [11,15] and b ∈ [0, 1]. Fig. 4(a) shows the SE (Spectral Entropy complexity) of system (1), with parameter a varying within [11,15] and b varying within [0.2, 0.6], SE complexity value is larger in all test range, the max value is 0.513. However, with parameter b varying within [0.6, 1], SE complexity value of system (1) falls instantaneously, even to 0. It means the dynamic behavior of the system change state from a chaotic state to a periodic state. The corresponding simulation result of C 0 complexity is shown in Fig. 4(b). The max value of C 0 complexity is 0.075, and its varying trend is basically the same to its varying trend of SE complexity. SE and C 0 complexity of the proposed system can be shown in Fig. 4(c) and 4(d). From the simulation result, when parameter a varying within [11,15] and b varying within [0.2, 0.6], the corresponding complexity spectrum become uniform and the complexity value has obviously gone up, the max value of SE complexity has already reached 0.581 and the max value of C 0 complexity is 0.093. The number of the periodic windows has been greatly reduced, and the complexity of the sequence has been greatly enhanced. By introducing two boosted control of cosine function, the periodic behavior has gradually changed into chaos when parameter b ∈ [0.6, 1], the primary lower complexity value has also change into the higher complexity value. All of these have a positive impact on its industrial application.

V. INFINITELY MANY COEXISTING ATTRACTORS
Different initial conditions can generate different chaotic orbits after many periods of evolution. However, it is difficult for attractors to overstep the constraints of the attractor basin and form large changes in shape or phase. For a chaos system of extremely multistability, the differential equations has an infinite number of solutions. The shift of the equilibrium point causes a shift in the phase of the attractor. Fixed system parameters of the proposed system, select the initial conditions is (x 0 ,1,z 0 ), x 0 and z 0 are variables. When a = 15, b = 0.4, set the initial values are (0,1,0), (0,1,π), (0,1,2π). . . ..(3π,1,3π), many coexisting attractors of periodic state are shown in Fig. 5(a). When a = 15, b = 0.1, and the initial values are (0,1,0). . . ..(3π,1,3π), many coexisting attractors of chaotic state are shown in Fig. 5(b). With two initial condition x 0 and z 0 varying simultaneously, the proposed system has extremely multistability. The phase position of attractors shifts periodically towards two directions of state variable x and z. It means this chaotic system can output chaotic signal with different amplitude at the same time. In comparison with system (1), the proposed chaos system can provide more ergodicity and flexibility for some engineering applications based on chaos theory.

VI. THE APPLICATION OF THE PROPOSED SYSTEM IN IMAGE ENCRYPTION A. THE DISCRETE SYSTEM
Due to the chaos degradation problem caused by the finite precision, the random signal of the continuous chaotic system cannot be directly applied to digital encryption. Therefore, in order to make the proposed chaotic system better suitable for digital encryption applications, it can often be further discretized. In which, by using Gaussian discrete algorithm, the proposed continuous chaotic system can be effectively discretized. The specific discrete systems equation are shown as follows: cos(2x(n))], y(n + 1) = y(n) + h[acos(2x(n))cos(2x(n))], where, a and b are control parameters, the step size h is set to 0.01 and n is number of cycles. Due to the characteristic of high pseudorandom, the chaos map is often used as a Pseudo Random Number Generator to generate the pseudo-random signal in engineering field.

B. CHAOTIC SECRET KEY
As shown in Fig. 6, a more than 300 bits-long secret key is applied in the design encryption algorithm. The design secret key is divided into 300 bits-long block S 1 and the secret key of initial value x(0), y(0), z(0). Where, we set secret keys x(0) = π, y(0) = 1, z(0) = π. Then S 1 is divided into S 11 , S 12 , S 13 , S 14 , S 15 , each of them can be divided into 60 bits-long. Especially in comparison with the single stable systems, the proposed new system of extremely multistability can generate infinitely many secret keys of initial value. This greatly increases the number of chaotic secret keys and thus increases the security of the encryption algorithm.

C. THE DESIGN OF ENCRYPTION ALGORITHM
Based on the pseudo-random values generated by the proposed discrete chaos system, the diffusion formula of adjacent pixels can be constructed. As Eq. (27) shown, each pixel value of the design diffusion formula depends on the value of two adjacent pixels and the proposed discrete chaos system. Then, doing the xor diffusion operation between each pixel value of the second design diffusion formula and the value generated by the proposed discrete chaos system. Then, doing the scrambling operation between three first pixels and three last pixels.
in which, p(n) is the pixel of the plaintext, x(n) is the value generated by the proposed discrete chaos system, E 1 (n) is the pixel value of the first diffusion algorithm, and E 2 (n) is the pixel value of xor diffusion algorithm, and a is the control parameter of the proposed discrete chaos system. To further improve the anti-deciphering ability of the algorithm, the second diffusion algorithm is designed, and it can be described as: where b is the control parameter, EE(n + 2) depends on two adjacent pixels, E 2 (n) and the proposed discrete chaos system, c 1 and c 2 are the parameter which depend on the values of interference parameters m 1 , m 2 , m 3 , m 4 . Then, do the third diffusion and scrambling operation: where EEE(n + 1) depends on the adjacent pixel EEE(n), the pixel of last diffusion EE(n) and the proposed discrete chaos system. Similarly, c 3 , c 4 depend on the parameter m 1 , m 2 , m 3 , m 4 and the value generated by the chaotic system. The images can be encrypted based on the design encryption formula, and the decryption algorithm is the inverse process of encryption algorithm.

D. THE ENCRYPTION AND DECRYPTION OF IMAGE
As Fig. 7 shown, based on the above proposed encryption algorithm, the plaintext image can be successfully encrypted and decrypted. In particular, the ciphertext image can effectively hide the core information of the original image. However, we cannot judge the performance of encryption only from the effect of ciphertext. The specific performance can be shown through the safety tests such as sensitivity of the secret key, histogram, information entropy and NPCR.

E. SENSITIVITY ANALYSIS OF THE SECRET KEY
Because the chaotic system is extremely sensitive to the initial state, extremely tiny variation of the initial secret key will generate a big difference after several iterations of the key sequence generator, and then the encryption effect can produce a big difference. It is because of this high sensitivity to initial values that chaotic systems are widely used in digital image encryption. In the design secret key, it is made up of 300 bits-long block S 1 and the secret key of initial values x(0), y(0), z(0) and then we will do the sensitivity tests on them. Leave the other parameters unchanged, Fig. 8(a)-(c) show two ciphered Hill image only by changing the secret key of initial value x(0) from π to π + 10 −14 . Fig. 8(d)-(f) show two ciphered Hill image only by changing the secret key of initial value y(0) from 1 to 1+10 −14 . Fig. 8(g)-(i) show two ciphered Hill image only by changing the secret key of initial value z(0) from π to π + 10 −14 . The experimental results show the good security of the encryption algorithm based on the design chaos secret key.

F. THE ANALYSIS OF HISTOGRAM AND CORRELATION
Histogram is widely used in many applications of the computer vision. It can detect changes in a video by marking significant edge and color statistical changes between frames. Image histogram is a histogram used to represent the brightness distribution in a digital image, it plots the number of pixels of each brightness value in the image. The histogram is a statistical collection of data and distributes the statistical results into a set of predefined bins. From the experimental results of Fig. 9, the histograms of the plain images contain the key information, and the histograms of the ciphered images using the design algorithm is fairly flat distributed, it can effectively hide the key information of the image. The correlation shows the strength between two adjacent pixels on the direction of horizontal, vertical and diagonal. It is one of most vital evaluation techniques in image encryption algorithm. The correlation coefficients test of Lena plain images and the ciphered images by using the proposed encryption algorithm are shown in Fig. 10. From the test result, the correlation of adjacent pixels in the plaintext image are higher, the concentration of the key information is high. However, the correlation of the ciphered image on the   direction of horizontal, vertical and diagonal are all relatively low, even they are relevant. And the specific data of correlation coefficients test are shown in Tab. 1.

G. ROBUSTNESS ANALYSIS
The ability to resist differential attack is often an important criterion to evaluate an algorithm. In general, tiny changes between adjacent pixels can cause major changes in image properties, NPCR (number of pixels change rate) and UACI (unified average changing intensity) are often used as a testing tool to measure this variation. The specific algorithm is described by: where k 1 , k 2 are two encrypted images generated by two plain images with a pixel difference and its size is L. Generally, the ideal value of NPCR is 0.9961 and the ideal value of UACI is 0.3346. The closer the values of NPCR and UACI obtained from the design encryption algorithm are to the ideal values, the higher the safety performance. That is to say, the proposed algorithm has have enough abilities to resist the differential attack. From the test results of NPCR and UACI in Tab. 2, the ability to resist differential attack of the design encryption algorithm is better. The information entropy of image represents a chaotic degree of pixels, it represents the randomness and uncertainty of image information. It is also an important evaluation standard of the safety performance of encryption algorithm. The    specific formula can be described as: in which, p(k i ) is the probability of k i . The ideal value of information entropy is 8, the closer the value of information entropy obtained from the design encryption algorithm are to the ideal value, the better the safety performance. The information entropy of the proposed encryption algorithm is shown in Tab. 3. From the comparison result, the pixels have higher chaotic degree and anti-deciphering ability.

VII. DSP IMPLEMENTATION
To make the chaotic signal better suitable for digital chaos applications, the proposed system based on the boosted control of cosine function can be implemented by Digital Signal Processing (DSP) technology. Select parameters a = 15, b = 0.1, and the starting condition (0,1,0). According to the implementation flow as shown in Fig. 11, the chaotic attractors on x-z plane can be digital implemented in Fig. 12. The experimental result of digital implement is basically the same to the analysis result of the numerical simulation, it shows the feasibility of theoretical analysis.

VIII. CONCLUSION
A new 3D chaotic system which can produce extremely multistability is reconstructed by the boosted control of cosine function. Base on its mathematical model, the Lyapunov exponents, Kaplan-Yorke dimension, bifurcation diagram and complexity can be simulated. In comparison with the original system, the proposed chaos system has better chaotic characteristic and higher complexity. Most of the cycle windows in the original chaos system became chaotic state, its chaos behavior was enhanced. Especially, the proposed system can generate infinitely many coexisting attractors due to the boosted control of cosine function. It means the proposed system can output more chaotic sequences of different amplitudes at the same time. These good performances have brought positive effects for its application in the industrial field. Next, we will continue to study the industrial applications based on this proposed chaotic system.