Enhancing Chaos Complexity of a Plasma Model through Power Input with Desirable Random Features

The present work introduces an analysis framework to comprehend the dynamics of a 3D plasma model, which has been proposed to describe the pellet injection in tokamaks. The analysis of the system reveals the existence of a complex transition from transient chaos to steady periodic behavior. Additionally, without adding any kind of forcing term or controllers, we demonstrate that the system can be changed to become a multi-stable model by injecting more power input. In this regard, we observe that increasing the power input can fluctuate the numerical solution of the system from coexisting symmetric chaotic attractors to the coexistence of infinitely many quasi-periodic attractors. Besides that, complexity analyses based on Sample entropy are conducted, and they show that boosting power input spreads the trajectory to occupy a larger range in the phase space, thus enhancing the time series to be more complex and random. Therefore, our analysis could be important to further understand the dynamics of such models, and it can demonstrate the possibility of applying this system for generating pseudorandom sequences.


Introduction
Further investigations in nonlinear dynamical systems have contributed to the development and understanding of numerous scientific phenomena. In particular, chaotic models as nonlinear systems have been derived from rules for describing complex behaviors. For instance, the Lorenz model is the 3D chaotic system proposed to represent the atmospheric convection [1], and the logistic map is the 1D discrete-time chaotic model proposed to describe the population growth [2].
Although chaotic systems are deterministic, which indicates that the whole future path of the system is uniquely governed by its initial condition, it is impossible to predict its long-term behavior [3,4]. These interesting systems have attracted considerable interest, especially in the area of telecommunication, security, laser, engineering, and many others [5][6][7][8][9]. Moreover, chaotic behavior can be observed in many disciplines of science such as physics, chemistry, biology, economy, and nature [10][11][12][13]. Particularly, it has been reported several plasma models exhibiting chaos. Some examples are as follows. The authors of [14] discovered the developed and weak chaotic motions in a magnetized dusty plasma model. The authors of [15] studied a new Lorenz-like system describing the interaction of three resonantly coupled waves in plasma. The authors of [16] presented the mechanical analysis and proposed a supremum bound of a plasma chaotic attractor.
The authors of [17] interpreted the dynamics of a mechanism plasma chaotic system. The authors of [18] proposed a method for finding hidden chaotic attractors in the plasma chaotic system.
The analysis of equilibrium points of chaotic systems has contributed to better understand their dynamics [19,20] and for determining their attractor type [21]. In general, there are two major types of chaotic attractors: the first type is called self-excited, and the second type is called hidden. The self-excited attractors can be found in chaotic systems that have a basin of attraction intersected with an unstable equilibrium [22]. Meanwhile, hidden attractors can be found in systems with a line of equilibria [23], curves of equilibria [24], no equilibria [25], or infinite number of equilibria [26]. The investigations of hidden and self-excited chaotic attractors have revealed unexpected results known as coexisting attractors or multistability, which has given the chaotic systems another dimension [27,28]. A gas laser model is the first system that shows multistability [29]; since that time many models with coexisting attractors have been presented [30][31][32][33][34]. From the application point of view, multistability provides large flexibility in the performance of nonlinear dynamical systems without changing system parameters [35,36]. Therefore, it can be used to induce switching between desirable attractors [27]. However, the multistability behavior can be harmful, especially in the design of a commercial device with particular properties where it is substantial to stabilize the required state in the existence of noise [37]. Additionally, the fluctuation of a chaotic system from chaos to the periodic state could break the chaos-based cryptosystems [30].
In this paper, we investigate the dynamics of a plasma perturbation, which consists of three coupled ordinary differential equations that contain three parameters [38]. The stability of equilibria reveals that the emergence of a chaotic attractor by the system is self-excited. Moreover, complex nonlinear behavior of transition from transient chaos to steady periodic can be observed by times series and phase space. To show the chaotic regions, which show no transient chaos, we use maximum Lyapunov exponents based on the contour plot. However, the system shows other complex behaviors of coexisting multiple symmetric attractors, which can be generated without adding forcing terms or controllers by increasing the power input. On this matter, boosting power input not only generates extreme multistability, but also can enhance the complexity and randomness of the system time series.
The remainder of this paper is organized as follows. Section 2 investigates the stability of the equilibria of the plasma perturbation model and determines its chaotic regions. Section 3 shows the ability of the plasma system for generating various multistability behaviors by increasing its power input. Section 4 illustrates the effect of increasing power input on the complexity of the system. In Section 5, we demonstrate the randomness of the system and show its high sensitivity to the initial conditions. Section 6 provides the conclusion of the paper.

The Mathematical Model
This paper concentrates on a low-dimensional model [38], that describes the dynamics of pellet injection in tokamaks. Mathematically, it is given by where p n is the normalized plasma pressure gradient, ξ n is the normalized magnetic field, and χ is the thermal diffusivity.
However, the corresponding 3D dynamical model can be found as follows, which was obtained by introducing the variables x 2 = χ anom χ 0 · ξ n and x 3 = p n . The constant parameters in the system (2) are defined as follows: d represents the dissipation/relaxation of perturbation, e represents the characteristic relation between the two heat diffusion coefficients, and h is the power input into the system.
We now give the basic features of the system (2) to interpret its complex behavior. First, the system (2) is symmetric about x 3 −axis, and this is due to that it is remained unaltered under the transformation (x 1 , . Second, it is dissipative only when the state space contraction, which is defined by (d + e + ex 2 2 ), is greater than zero. Finally, three different equilibrium points-P 1 (0, 0, h), P 2 (0, √ h − 1, 1), and P 3 (0, − √ h − 1, 1)-can be obtained by considering the left hand sides of the three equations in the system (2) equal to zero.
For studying the stability of the system (2) at the equilibrium point P 1 , we first obtain the corresponding eigenvalues, which are given by A continuous dynamical system is unstable when its eigenvalues have positive real parts. Obviously, λ 1 is positive only when the parameter e is less than zero. Meanwhile, λ 2,3 have positive real parts if h > 1. Moreover, the characteristic equation at P 2,3 is given by Suppose that the parameters d, e, h are greater than zero. Then, the eigenvalues at P 2,3 have negative real parts when Thus, the equilibrium points P 2,3 are unstable when one of the above inequalities does not hold. In other words, the system (2) at P 2,3 is unstable if and only if

Self-Excited Chaotic Attractor and Complex Transient Chaos
When the parameters values set as d = 0.5, e = 0.1, h = 2, Figure 1a-c depicts the time series, largest Lyapunov Exponents, and the phase space of the system (2), respectively. For the numerical simulation of the system (2), we use Matlab 2019a program to plot the time series, the phase spaces, and largest Lyapunov Exponents, in which the standard Runge-Kutta fourth-order scheme was utilized via Matlab to obtain the system solution and then plot the time series and the phase spaces, whereas Wolf algorithm [39] was utilized to calculate the largest Lyapunov Exponent in Matlab. However, it can be seen that the system exhibits chaotic attractor for a long duration time when the starting point is from a neighborhood of the equilibrium point P 1 . As all equilibria of the system (2), including P 1 , are unstable for d = 0.5, e = 0.1, h = 2, the generated chaotic attractor in Figure 1 is self-excited. On the other hand, transient chaos usually appears when the initial conditions are chosen near the boundary of the basin of attraction or when the system is very close to a bifurcation point. Thus, by choosing the parameters values as d = 0.5, e = 0.1, h = 1.179, which make the system is close enough to a bifurcation point, we observe that the system (2) exhibits a complex chaotic transient, as shown in Figure 2. In this figure, we have plotted the time series and the corresponding phase portraits for a specific set of the system parameters and with two different sets of initial conditions, which have been chosen near the boundary of the basin of attraction. As can be seen, for each set of the initial conditions, all variables of the system (x 1 , x 2 , x 3 ) show a chaotic behavior in short duration time, and subsequently those variables show asymptotic orbit for the rest of the time.

Chaotic Regions of the Plasma Model
Chaos is usually identified by one of the fundamental algorithms, which is Lyapunov exponent. The Lyapunov exponent is a quantitative measure employed to estimate the divergence or convergence of adjacent trajectories in a system. Here, we calculate the largest Lyapunov exponent when two parameters are varying simultaneously, which can provide a wider vision of the chaotic regions of the system. Therefore, three different sets of system parameters are selected to depict the largest Lyapunov exponent with a step size of 0.02, as shown in Figure 3. Based on this figure, it can be observed that the system (2) exhibits three different behaviors, which can be determined as follows. • Figure 3a illustrates that the chaotic behavior appears in the brown, yellow, and green colors regions. Meanwhile, the quasi-periodic and periodic behaviors appear in cyan and blue colors regions, respectively. • Figure 3b,c shows that the chaotic behavior of the system (2) appears in the brown color region, and the quasi-periodic behavior appears in the green color region. Meanwhile, the cyan and blue colors regions indicate to the periodic behavior of the system.

The Plasma Model with Coexisting Symmetric Attractors
In this section, we study the dynamical behaviors of the plasma model (2) by injection more power input to the system. This can be done by adding a constant parameter A to the third equation and fixing the diffusion parameters e. As a result, the modified system is defined as follows, where A is a positive parameter. Obviously, adding A to the third equation of the system (2) can keep the symmetry of the system about x 3 −axis intact. Moreover, the system (9) has the same number of equilibria as in the system (2). These equilibria are given by The stability of P 1 , P 2 , P 3 is determined by the following proposition.

Proposition 1.
The following statements hold if the parameters d, h, e, A are greater than zero.

The Coexistence of a Symmetric Pair of Attractors
A chaotic system, which exhibits different attractors for a particular set of system parameters, is considered as highly sensitive to its initial conditions. Therefore, this section discover the presence of coexisting attractors after the injection of more power input to the system.
By selecting A as control parameter for the range 0.19 ≤ A ≤ 0.3, the largest Lyapunov exponents and multistable bifurcation diagram of the system (9) are, respectively, plotted in Figure 4a,b when the step size is equal to 0.0001. By observing such the two figures, the system is simulated with two different initial conditions simultaneously, in which the blue and red orbits begin with the initial conditions (1, 1, −1) and (2, 1, −1), respectively. Obviously It is crucial here to plot the basins of attraction to provide foresight about the behavior of the system with a wide range of initial conditions. Therefore, the basins of attraction of the system in two different scenarios are plotted: the parameters region where two chaotic attractors coexist, and the region of the parameter where two quasi-periodic at-tractors coexist, as illustrated in Figure 5a,b, respectively. For both scenarios, the basins of attractions agree with the bifurcation diagram in which the system exhibits two symmetric coexisting attractors. This can be further illustrated by depicting the corresponding phase portraits, as shown in Figure 6. It can be seen that the system with A = 0.22 shows the coexistence of symmetric wings of strange attractors in which these attractors have formed a shape that resembles a butterfly, as shown in Figure 6a,b. Besides that, the coexistence of symmetric wings of double-period periodic attractors can be observed with A = 0.24, as shown in Figure 6c,d. However, Table 1 lists the Kaplan-Yorke dimension (DKY) and Lyapunov exponents (LE) of the attractors that are shown in Figure 6. Here, we use the Wolf algorithm [39] to calculate the Lyapunov exponents with t = 10 4 .

The Coexistence of Many Symmetric Quasi-Periodic Attractors
To examine the effect of increasing the power input for the system, Figure 7 depicts the coexisting bifurcation diagram with respect to parameter A. In this figure, we set the parameters d = 0.5, e = 0.1, h = 1.95 and 0.39 ≤ A ≤ 0.43. Meanwhile, the system is simulated with several sets of initial values, which are selected as (K, 1, −1). When K is set to −6, −3, −1, 0, 2, 4, 5, the coexistence of seven quasi-periodic attractors is observed mainly within the range A ∈ [0.39, 0.41]. This interesting nonlinear phenomenon can be further illustrated by plotting different projections of the system phase portraits, as shown in Figure 8.

Sample Entropy Algorithm
This section studies the complexity performance of systems (2) and (9) by one of fundamental algorithms, which is Sample Entropy (SamEn) [40]. SamEn is derived from approximate entropy by Richman et al. to estimate how much extra information is required to predict the (t + 1)th output of a trajectory using its previous (t) outputs. Larger SamEn values reflect a lower degree of regularity, which means that higher complexity performance.

1.
Reconstruction: the time series can be reconstructed as follows where m is embedding dimension, and τ is time delay.

2.
Counting the vector pairs: For a given tolerance parameter r, let B i be the number of vectors X j such that here, d[X i , X j ] is the distance between X i and X j , which is defined as 3.
Calculating θ m (r): According to the obtained number of vector pairs, we can get then calculate θ m (r) by 4. Calculating SamEn: Repeating the above steps we can get θ m+1 (r), then SamEn is given by For m = 2 and r = 0.2× Standard Deviation, we depict the complexity analysis results of the system (9) using Matlab 2019a program with two different sets of the system parameters, as shown in Figure 9a,b. There are two parameters simultaneously varying with the step size equal to 0.01, in which one of them is always A. Choosing A as a varying parameter can demonstrate the effect of increasing the power input on the complexity of the system. As can be seen in Figure 9, the complexity of the system is enhanced with the increase of A. In other words, boosting power input to the system can decrease its degree of regularity.
To perceive the effect of decreasing the degree of regularity on the dynamics of the system, Figure 10 plots the phase portraits comparisons with two values of the parameter A. In this figure, the trajectories with lower complexity (red color) are generated by the system (9) for A = 0. In other words, these trajectories are the numerical solutions of the system (2), whereas the trajectories with higher complexity (blue color) are generated by system (9) for A = 0.105. Thus, it can be observed that enhancing the complexity of the system (2) could change its nature from order to chaos or even expand its chaotic oscillations to occupy a wider region in the phase space, as shown in Figure 10a

Performance Evaluations
This section evaluates the sensitivity and randomness of the chaotic sequences generated by the system (9) by cross-correlation coefficient and NIST-800-22 test.

Cross-Correlation Coefficient
As can be seen in Section 3, in the multistability regions, system (9) is sensitive to the initial conditions. However, it is important to examine the sensitivity performance of the system in the single stability regions, which are more desirable in cryptography applications [30]. Therefore, this subsection deals with the two cases that have appeared in Figure 10. The sensitivity of system (2) and (9) can be estimated by the cross-correlation coefficient (CCF), which is given by where A(α) represents the mean value of the time series α t , meanwhile A(β) represents the mean value of the time series β t . When the CCF(α t , β t ) result is close to 0, then it can be indicated that these two time series are diverging. When d = 0.5, e = 0.1, h = 2, the sensitivity of systems (2) and (9) is illustrated in Figure 11a,b, respectively. In these two figures, the CCF is calculated between the original time series (S 1 ), and the modified one (S 2 ) which is created by contaminating one of the initial conditions of the system via a tiny error E = 5 × 10 −5 . Moreover, the sensitivity of the systems (2), and (9) are respectively depicted in Figure 11c,d, where the parameters are set as d = 1, e = 0.1, h = 2, and A = 0.105. The system (9) has a higher sensitivity to its initial conditions.

Chaos-Based Cryptographic Pseudo-Random Number Generator (PRNG)
A chaotic system with complex nonlinear behaviors is an appropriate source for producing PRNG. Generally, the randomness of the generated PRNG by chaotic systems is dependent on the characteristic properties of these systems such as the sensitivity of initial state, ergodicity, and unpredictability. As system (9) is extremely sensitive to its initial conditions and has high complexity performance, this system could be suitable for chaos-based PRNG.
To further examine the ability of the system (9) for generating valid PRNG, we propose a simple strategy, which directly uses the chaotic sequences as pseudorandom numbers. For the state variables of system (9), {X(i), Y(i), Z(i)|i = 1, 2, 3 . . . }, we convert each of their values to 32-bit binary stream using IEEE 754 float standard, as shown in Figure 12. Subsequently, the digital numbers in each binary stream are extracted from 22nd to 32nd (for X(i), and Y(i)), and from 23nd to 32nd (for Z(i)) to compose PRNG of 32-bits.
The validity of the generated PRNG by the system (9) can be checked using NIST-800-22 [41]. In NIST-800-22, there are 16 different statistical tests. However, a random sequence can pass the test if the corresponding p-value is greater than the experimental significance level of α. In our experiment, we generate a binary sequence of 10 8 bits, which means that the obtained p-value should be greater than α −1 in each test. Figure 13 illustrates the NIST SP800-22 test results of PRNG that generates by the system (9). As can be seen, all the p-values are greater than α −1 , which demonstrates the high randomness of the PRNGs using the system (9).

Conclusions
The dynamics of the 3D plasma perturbations model, which was proposed by Constantinescu et al. to describe the pellet injection, was investigated by time series, Lyapunov exponents, bifurcation diagrams, and basins of attraction. Besides that, the stability of its equilibria has been analyzed. The mathematical analysis has shown that the system has three equilibria in which there is always an unstable equilibrium point, which indicates that the attractors type of this system is self-excited. To discover the chaotic regions that show no transient chaos, the two-dimensional MLE has been applied. Furthermore, we have conducted comprehensive research on the dynamics, complexity, and randomness of the system after boosting power input. Simulation results have demonstrated that the system could produce multiple coexisting attractors, high complexity, and randomness.

Conflicts of Interest:
The authors declare no conflicts of interest.