**2. Chaotic Flow with a Hyperbolic Sinusoidal Function**

In the search for chaos flows with hyperbolic sinusoidal function, we study the form of a three-dimensional chaotic structure as:

$$\begin{aligned} \dot{x} &= a\_1 x + a\_2 y + a\_3 z + a\_4 x y + a\_5 x z + a\_6 y z + b, \\ \dot{y} &= a\_7 x y + a\_8 x z + a\_9 y z + a\_{10} \sinh(y) + c, \\ \dot{z} &= a\_{11} x + a\_{12} y + a\_{13} z + a\_{14} x y + a\_{15} x z + a\_{16} y z \end{aligned} \tag{1}$$

where *x*, *y* and *z* denote the system states; *a*1, ... , *a*<sup>16</sup> indicate the coefficients of the terms; *b* and *c* are two scalars which define the chaos behavior.

A computer examination is executed investigating millions of combinations of different forms, various initial states and different constant values, looking for dissipative cases for which the largest Lyapunov exponent is bigger than 0.001. The system is in chaos state if the largest Lyapunov exponent is bigger than zero, and the system is in steady period state if the largest Lyapunov exponent is smaller

than zero [38,39]. Therefore, in the present work, a three-dimensional chaotic flow is reported which is specified by: .

$$\begin{cases} \dot{x} = x - ayz + b, \\ \dot{y} = xz - \sinh(y) + c, \\ \dot{z} = x \end{cases} \tag{2}$$

where the parameter *a*<sup>6</sup> in (1) has been denoted as parameter *a* in system (2).

Next, three different scenarios depending on the values of system's (2) parameters *b*, *c* are discussed in details.

## *2.1. Scenario A: Line of Equilibria*

If *b* = *c* = 0, the chaotic flow (2) will have a line of equilibria, i.e., *EA* = [0, 0, *z*∗] *<sup>T</sup>*, where *z*<sup>∗</sup> is the equilibrium point value in *z* axis and *T* means the transpose of the vector. To analyze the state trajectories in the vicinity of the equilibrium point, the Jacobian matrix is obtained from (2) as:

$$J = \begin{bmatrix} 1 & -az & -ay \\ z & -\cosh(y) & x \\ 1 & 0 & 0 \end{bmatrix} \tag{3}$$

For equilibrium point *EA* = [0, 0, *z*∗] *<sup>T</sup>*, the Jacobian matrix is defined as

$$J = \begin{bmatrix} 1 & -az^\* & 0 \\ z^\* & -1 & 0 \\ 1 & 0 & 0 \end{bmatrix}. \tag{4}$$

Therefore, the eigenvalues of the linearized system are achieved as

$$|\lambda I - f| = \begin{vmatrix} \lambda - 1 & az^\* & 0 \\ -z^\* & \lambda + 1 & 0 \\ -1 & 0 & \lambda \end{vmatrix} = \lambda \begin{pmatrix} \lambda^2 + az^{\*2} - 1 \end{pmatrix} = 0 \Rightarrow \lambda\_1 = 0, \lambda\_{2,3} = \pm \sqrt{1 - az^{\*2}}.\tag{5}$$

The equilibrium is a saddle node for <sup>−</sup> <sup>1</sup> √ *<sup>a</sup>* <sup>&</sup>lt; *<sup>z</sup>*<sup>∗</sup> <sup>&</sup>lt; <sup>1</sup> √ *a* . The equilibrium point is an unstable node for *<sup>z</sup>*<sup>∗</sup> = <sup>±</sup> <sup>1</sup> √ *a* . For the values *z*<sup>∗</sup> > <sup>1</sup> √ *<sup>a</sup>* and *<sup>z</sup>*<sup>∗</sup> <sup>&</sup>lt; <sup>−</sup> <sup>1</sup> √ *a* , since one eigenvalue is zero and two eigenvalues are imaginary, the stability of the equilibrium point cannot be determined by this method; the equilibrium point may be stable, unstable or marginally stable.

If one design parameter is varied and the norm of the state variables vector is plotted for finding the fixed points of the system versus the changing parameter, finally the bifurcation diagram is obtained [40]. In the bifurcation diagrams, the fixed points maybe disappear, appear, or change their stability nature when the design parameter is changed. Those variations may occur even for infinitesimal changes in the parameter. Bifurcation diagram is used for the stability analysis of a dynamical system [41,42]. Moreover, the Lyapunov exponents spectrum makes it possible to qualitatively quantify a local property with respect to the attractor's stability. The positive/negative values of the Lyapunov exponents can be observed as a measure of the averaged exponential divergence/convergence of neighborhood trajectories [43,44]. The bifurcation diagram for *y*, when the states cut the plane *z* = 0 with d*z*/d*t* < 0*,* as well as the spectrum of system's Lyapunov exponents (*LEi*, *i* = 1, 2, 3), by varying the value of *a* to explore the dynamical form of system (2), while keeping the initial states as [*x*0, *y*0, *z*0]=[2, 0.2, 1], are depicted in Figure 1. Therefore, the suggested structure (2) is integrated via the classical Runge-Kutta integration algorithm [45], numerically. For all of the parameters, the simulation calculations are executed via the parameters and variables in extended precision mode. In addition, the spectrum of the Lyapunov exponents are found via the Wolf's algorithm [46].

(**b**)

**Figure 1.** (**a**) Bifurcation diagram, (**b**) Lyapunov exponents spectrum of dynamics (2), when changing *a* from 1.85 to 2.4, and *b* = *c* = 0.

The dynamics (2) shows a chaotic attractor, for *a* = 2 (Figure 2), and a limit cycle of Period-1, for *a* = 2.35 (Figure 3). The spectrum of Lyanpunov exponents (Figure 1b) approves the dynamic behavior of the system as it has been revealed via bifurcation diagram.

## *2.2. Scenario B: No Equilibrium Point*

If *b* - 0, *<sup>c</sup>* <sup>=</sup> 0 and by keeping *<sup>a</sup>* <sup>=</sup> 2, for obtaining the equilibrium point, we solve . *<sup>x</sup>* <sup>=</sup> 0, . *y* = 0 and . *z* = 0, that is

$$\begin{cases} x - 2yz + b = 0, \\ xz - \sinh(y) = 0, \\ x = 0. \end{cases} \tag{6}$$

Consequently, the chaotic flow has no equilibrium point in this case. Therefore, it belongs to the category of chaotic systems containing hidden attractors.

Taking the bifurcation diagram of *y* (Figure 4a), along with the Lyapunov exponents spectrum (Figure 4b) by changing *b* for 0 < *b* < 0.005, in order to explore the dynamics (2), for initial conditions [*x*0, *y*0, *z*0]=[2, 0.2, 1], interesting dynamical behavior has been investigated. As it is obtained from bifurcation diagram (Figure 4a), the system passes from a chaotic region, for *b* ∈ [0, 0.075), to a periodic one as the parameter *b* increases.

**Figure 2.** Strange chaos attractor for *a* = 2 and *b* = *c* = 0 in Scenario A.

**Figure 3.** Limit cycle of period-1 for *b* = *c* = 0 and *a* = 2.35 in Scenario A.

**Figure 4.** (**a**) Bifurcation diagram, (**b**) spectrum of Lyapunov exponents of (2), when changing *b* from 0 to 0.01, for *a* = **2**, *c* = **0**.

The strange attractors of the system (2) are displayed for *a* = 2, *b* = 0.005 and *c* = 0 in Figure 5. In this case, the Lyapunov exponents are *LE*<sup>1</sup> = 0.14107, *LE*<sup>2</sup> = 0, *LE*<sup>3</sup> = −1.33835, which confirmed the chaotic behavior of the system (2). The Kaplan-York dimension of the chaotic flow is *DKY* = 2.1054. Besides, the Poincaré map in *x*-*y* plane presents the folding properties of chaos when *z* = 0 with d*z*/d*t* < 0 (Figure 6).

**Figure 5.** Strange chaotic attractors for *a* = 2, *b* = 0.005 and *c* = 0.

**Figure 6.** Poincaré map of chaotic system (2) in the *x–y* plane, for *a* = 2, *b* = 0.005 and *c* = 0.

#### *2.3. Scenario C: Self-Excited Attractor*

If *b* = 0, *c* - 0 and *<sup>a</sup>* <sup>=</sup> 2 this chaotic flow has only one equilibrium *EC* <sup>=</sup> 0,*sinh*−1(*b*), 0*<sup>T</sup>* . For the equilibrium point *EC*, the Jacobian matrix is found as:

$$J = \begin{bmatrix} 1 & 0 & -2\sinh^{-1}(b) \\ 0 & -\cosh(\sinh^{-1}(b)) & 0 \\ 1 & 0 & 0 \end{bmatrix} . \tag{7}$$

Then, the eigenvalues of the linearized chaotic flow are obtained as:

$$\begin{aligned} |\lambda I - f| &= \begin{vmatrix} \lambda - 1 & 0 & 2\sinh^{-1}(c) \\ 0 & \lambda + \cosh(\sinh^{-1}(c)) & 0 \\ -1 & 0 & \lambda \end{vmatrix} \\ &= \left(\lambda + \cosh(\sinh^{-1}(c))\right)\left(\lambda^2 - \lambda + 2\sinh^{-1}(c)\right) = 0 \\ &\Rightarrow \lambda\_1 = -\cos(\sinh^{-1}(c)) = -\sqrt{1 + c^2}, \lambda\_{2,3} = \frac{1 \pm \sqrt{1 - 8\sinh^{-1}(c)}}{2}. \end{aligned} \tag{8}$$

For *c* > 0.1253, the eigenvalues of the chaotic flow are found as λ = − √ 1 + *c*2, <sup>1</sup>±*i*<sup>ω</sup> <sup>2</sup> , and the equilibrium point is a saddle focus. For *c* < 0.1253, the eigenvalues of the chaotic flow are obtained as λ = − √ 1 + *c*2, <sup>1</sup><sup>±</sup> √ Δ <sup>2</sup> , and the equilibrium point is a saddle node.

Figure 7 depicts the bifurcation diagram of variable *y* as well as the spectrum of Lyapunov exponents by varying *c,* for 0<*c*<0.05, to explore the dynamics (2), for initial states [*x*0, *y*0, *z*0]=[2, 0.2, 1]. It is shown from bifurcation diagram (Figure 7a) that the system passes from a chaotic region for *c* ∈ [0, 0.0285) to a periodic one as the parameter *c* increases. The respective spectrum of Lyapunov exponents to parameter *c* displays the aforementioned system's (2) dynamical behavior for *a* = 2 and *b* = 0.

The strange attractors of (2) for *a* = 2, *b* = 0, and *c* = 0.02 are demonstrated in Figure 8. For these parameter's values, the Lyapunov exponents are *LE*<sup>1</sup> = 0.09822, *LE*<sup>2</sup> = 0, *LE*<sup>3</sup> = −1.27669, which confirmed the chaotic behavior of system (2). The Kaplan-York dimension is *DKY* = 2.07699. Furthermore, the Poincaré maps in *x*–*y* plane, when *z* = 0 with d*z*/d*t* < 0 (Figure 9) presents the folding features of chaos.

**Figure 7.** (**a**) Bifurcation diagram, (**b**) Lyapunov exponents spectrum of (2) when changing *c* from 0 to 0.05, for *a* = **2**, *b* = **0**.

**Figure 8.** Strange attractors for *a* = 2, *b* = 0 and *c* = 0.02.

**Figure 9.** Poincaré map of system (2) in the *x-y* plane, for *a* = 2, *b* = 0 and *c* = 0.02.
