Multistability Route in a PWL Multi-scroll System through Fractional-Order Derivatives

Recently, it was discovered that the use of fractional derivatives induces the occurrence of multistable states in PWL systems with multiple scrolls. In this paper, we show the path of the system from monostable to multistable behavior by reducing the integration order. Using bifurcation diagrams, which describe the evolution of the attractor under the modiﬁcation of the bifurcation parameter, and basins of attraction, we show that the system has a path to multistability, revealing the mechanism for the emergence of the behavior. The resulting dynamics shows the coexistence of up to n + 2 attractors, where n is the number of scrolls generated by the integer-order system. Moreover, the existence and uniqueness of the dynamics is proved by Poincar´e sections showing that the behaviors exist and coexist and are not a section of the same solution.


Introduction
The coexistence of possible final states for a certain fixed set of parameters, where the final state depends on the system's initial conditions, is known as multistability [1]. This phenomenon has been a current topic in scientific development due to two factors; the high frequency of appearance in natural systems and the complexity involved in controlling such behavior to take advantage of its properties in technological applications [2][3][4][5]. Multistability is a phenomenon mainly associated with dissipative dynamical systems, having as main mechanisms for its appearance the coupling of systems, state feedback, systems with delays, the appearance of multistability due to parametric forcing, induced by noise, Gavrilov-Shilnikov -Newhouse sinks, and the appearance of coexisting states in slightly dissipative systems [1]. In a PWL system, linear and nonlinear components interact, leading to the special properties of generating multi-scroll attractors [6][7][8][9][10]. Various approaches to designing multi-scroll systems include adding breakpoints to Chua's function, using hysteresis functions, applying periodic functions, by using Piece-Wise Linear functions (PWL), with results ranging from 1-D to 10-D, and generating self-excited and hidden attractors.
It is possible to find the appearance of coexisting states in PWL systems [11][12][13][14][15][16][17][18], for example, through the coupling of oscillators with different number of scrolls, giving rise to both a generalized synchronization and an itinerant state capable of living in the entire space of the system with a greater scroll number, but generating the minimum amount induced in the coupling of systems [19,20]. At the same time, this type of system presents multistable states through the parametric modification that generates attractors of a single-scroll [11][12][13], or through the construction of ad-hoc nonlinear functions that induce the coexistence of possible final states [21], considering the described methods from the perspective of traditional computation (integer-order calculus).
With the recent increase of studies and experiments with fractional order systems, the possibilities of finding new behaviors and better descriptions of natural phenomena are a recurring theme in the literature [22][23][24][25][26][27][28][29][30]. However, the use of this numerical tool has been neglected because it is used as a dynamical validation mechanism and the effects and physical implications associated with the use of fractional order derivatives are ignored. In [31], the authors provide a review of several papers highlighting various physical phenomena where the use of fractional order derivatives allows us to better describe the mechanisms by which the estimated behavior is generated by traditional derivatives, and emphasize the relevance of using fractional order techniques with nonlinear oscillators.
Under the premise of better understanding the effects and physical mechanisms in fractionalorder systems, it was recently published that the use of fractional-order derivatives in PWL systems induces the occurrence of coexisting states that breakup the basin of attraction of a monostable system into n possible basins, where n is the number of scrolls associated with the integer-order system [23]. In such work, the phenomenon has only been reported and it has been restricted to the study of fractional derivative values where such behavior is possible. Further explanations for the mechanism that triggers the multistable behavior have not been explored. The relevance of the study lies both in the novelty of the discovered behavior and in the application of the jerk-inspired system to various physical phenomena such as Brownian motion [32,33].
Considering that the occurrence of multistability in a PWL system has been reported, but no further analysis or clarification of the physical process involved in the appearance of the behavior, this paper analyzes the mechanism by which a multiscroll system fragments its basin of attraction due to the fractional calculus approach. As a result, a path of the system from monostability to multistability is discovered when a nonlinear jerk oscillator with fractional order derivatives is studied. The obtained results show that the system exhibits smooth transitions when the fractional order is reduced. Here, the system first divides its attraction region into n + 2 sections for a small number of values when the integration order is reduced, which opens the possibility of finding the same number of coexisting states, and later converges to only n possible attractors before the dynamics is finally canceled due to system stability.
The rest of the work is as follows: The necessary preliminaries to describe the oscillator and the fractional-order derivative used are discussed in Section 2. The multistability route through fractional-order characterized by bifurcation diagrams and basins of attraction is addressed in Section 3. Section 4 describe the Poincaré analysis and its simplications. The discussion and extension of the results to systems with greater number of scrolls is presented in Section 5. Conclusions are devoted to the last section of the paper.

Fractional-Order Calculus
With the advent of fractional order theory, developed nearly 300 years ago, the study of chaotic systems through fractional order derivatives has been extensively studied by the scientific community [34][35][36], due to the memory and inheritance properties presented by several natural materials and processes [37,38].
Derivatives and integrals of integer order are generalizations of the analogous derivatives of fractional order. Various definitions of fractional order derivatives can be found in the literature, the most common being Riemann-Liouville and Caputo operators [39,40]. The Caputo definition of fractional derivative is described by: with n = ⌈q⌉, where Γ represents the Gamma function defined as follows: A commensurate fractional-order timeinvariant system can be described in a general form as in Eq. (3): which subject to initial conditions, turns into where n 1 , n 2 , . . . , n k are rational numbers, such that n k > n k−1 > · · · > n 1 > 0, n j − n j−1 ≤ 1 for all j = 2, 3, . . . , k and 0 < n 1 ≤ 1. The least common multiple of the denominator of n 1 , n 2 , . . . , n k is defined by M , and set q = 1/M and N = M n k . Then Eq. (3) can be expressed as in [39]: considering initial conditions as: Furthermore, this linear time-invariant system can be expressed in a matrix form as in Eq. (7).
where x ∈ R n is the state vector, A ∈ R n×n is a linear operator, and q is the fractional commensurate derivative order 0 < q < 1. Unlike integer-order systems, the stability of fractional-order systems depends on the derivative order (q), which generates a stability region, as shown in Figure 1 [41]. It is worth noting that the stability of an equilibrium point can be modified by the fractional order. Considering a general jdimensional system of fractional order such as the one in Eq. (7), with λ j -eigenvalues, it is possible to define the stability of the system by its eigenvalue analysis, and classify it as follows [39]: ❼ The system is stable, if and only if |argλ j | ≥ qπ 2 . ❼ The system is asymptotically stable, if and only if |argλ j | > qπ 2 . ❼ The system is unstable, if and only if |argλ j | < qπ 2 , for at least one eigenvalue. The interest of this work is to preserve unstable equilibrium points, hence the system in Eq. (7) must have at least one eigenvalue in the unstable region, i.e., if the system does not satisfy the locally asymptotic stable condition:

PWL Multi-scroll System
Consider a PWL multi-scroll system based on the jerk equation [12,[42][43][44][45], which satisfies all conditions to be considered as UDS I for the control parameter space α ∈ (0 − 1), since it has a real positive eigenvalue, a pair of complex conjugate eigenvalues with negative real component, where the eigenvalue sum function is negative.
where D q is the fractional order Caputo operator. The control parameter α is a positive constant defining the eigenvalues of the system, where α ∈ R. It is possible to generate an attractor with multiple scrolls by using a commutation law f (x), in this case a Piece-Wise Linear (PWL) function, as shown in Eq. (10) [23,44], whose purpose is to control the visitation in the different domains containing the equilibrium points of the system, which is achieved by the coexistence of a large number of unstable one-spiral trajectories.
As mentioned in [22,23], the system has bounds on the value of the fractional derivative that can be modeled without the system converging to a stable equilibrium point. Considering a fixed value at α = 0.45, the minimum order of integration before stabilizing the dynamics of the system is q c = 0.8958. Figure 2 shows both the integer order attractor of Eq. (9-10) and the bifurcation diagram of the oscillator when the integration order is changed by considering random initial conditions. These results are identical to those shown in [23] computed using the Adams-Bashforth-Moulton (ABM) method for fractional order systems [46]. Bifurcation diagram, Figure  2(b), shows the multistable region previously identified by our research group, where it is possible to find the coexistence of n single-scroll attractors, being in this case n = 3.

Multistability Route through Fractional-Order
Considering the results shown in the bifurcation diagram of Figure 2(b), where only the occurrence of multistability in the multi-scroll system is indicated due to the decrease of the integration order, the transition of the system between the multistability and monostability is investigated. Using the same integrative method as in the previous sections, the range of integration order studied is increased from q ∈ [0.8958, 0.8978], multistable reported region, to q ∈ [0.896, 0.899]. By monitoring the dynamic changes of the attractor, taking the integration order as bifurcation parameter and a fixed condition, the steady-state dynamics achieved by the system for these parameters is captured, then the last value of the time series is taken and injected as the new initial conditions of the model for a small change in the bifurcation parameter [47,48]. For this purpose, the bifurcation parameter is increased to control the multistable states at the beginning. The process is repeated until the upper integration limit is reached. This method of creating the bifurcation diagram, in contrast to the consideration of random initial conditions, that reveal the global dynamics of the system, allows analyzing the evolution of the attractor under the modification of the bifurcation parameter and provides the opportunity to evaluate the transitions in the limit to monostability. As result, the behavior depicted in the bifurcation diagram in Figure 3 Figure 3 shows the bifurcation diagrams obtained by analyzing the dynamic evolution of the attractor under the modification of the integration order, considering an integration step size ∆ q = 6.39e −6 . The individual processing of the time series allowed the identification of two new attractors, the coexistence of two double-scroll dynamics once it is "inside" of the reported multistability region (q ∈ [0.8958, 0.8978]), together with the single-scroll multistability response. These states were not reported in previous studies for both integer-order and fractional-order PWL systems due to the method used to construct the bifurcation diagram.
Note that each of the bifurcation diagrams shown in Figure 3 is generated under a different initial condition and that a small change in these conditions can lead to a different state in one of the multistable regions of the system. By analyzing the evolution of the attractor as the derivative order is gradually changed, it is possible to see the path of multistability as the system moves from a single-scroll attractor to visit two domains of the system and finally converge to the monostable dynamics with 3 scrolls. Moreover, the colors associated with the branches of each diagram refer to the attractor obtained for these orders of integration. Figure 4 shows the attractors included in the system by changing the integration order. Figure 4 shows all possible attractors resulting from the bifurcation diagrams of Figure 3. To validate the behaviors shown, the Maximum Lyapunov Exponent (MLE) was calculated from the time series as defined in [49], resulting in the following values:

Poincaré Analysis
Since the phenomenon presented here is limited to numerical proofs, and in order to guarantee that the above described multistable solutions coexist in the space of states and not that one is a complement to another (i.e., that the multistable state of 1-scroll attractor is part of the solution of the double-scroll attractor, which has not been properly quantified), Poincaré planes are used to graphically represent the crossing points between the domains where the scrolls are generated by the switching surfaces [45]. Consider the PWL system described in Eq. (9), where f (x) is given by b i ∈ R 3 , i = 1, 2, . . . , n constant vectors inducing n polytopic partitions of the state space, and let Ω 1 , ..., Ω n be the associated domains to the PWL function such that n i=1 Ω i = R 3 and Ω i (Ω j ) 0 = ∅, where the notation (Ω j ) 0 denotes the interior of the domain Ω j . Moreover, in each domain Ω i ⊂ R 3 the system has an equilibrium point χ * i = −A −1 b i and evidently there exist as many equilibria as domains Ω i . Under this assumption, it is possible to implement a Poincaré plane at the commutation generated by b i ∈ f (x). Under this assumption, the Poincaré plane supervises the trajectory transitions from Ω i → Ω j or the return of the trajectory to the same domain Ω j → Ω i .
To illustrate the construction of the Poincaré section and the information collected, consider the multistable attractor located around the origin, as shown in Figure 4  with the switching surface from Ω i → Ω j , and Ω j → Ω i , respectively. The color in the Poincaré sections is used only here to distinguish the points associated with each commutation. Figure 5(a,b) illustrates the 1-scroll multistable attractor, the location of the Poincaré section, such as the recurrence points. Figure 5(c,d) shows the location of the intersection points obtained for the attractor shown in Figure 5(a,b). For this case, is possible to see that the crossing points have a kind of symmetry, that is, the exit points of the trajectory from the central domain are represented by black crosses and red circles, and the entry points into the central area are represented by black circles and red crosses.
Since the goal of analyzing multistable behavior with a Poincaré section is to identify the uniqueness of the solution and analyze the behavior of each attractor at the switching surfaces, the crossing points in each plane of all multistable attractors are determined and the results are shown in Figure 6. The color in each figure corresponds to the color of the attractor. Figure  6(a) shows the intersections of the trajectories at commutation in S 1 , while Figure 6(b) shows the transitions of the trajectories in S 2 . Note that each plane is crossed by only three attractors, for plane S 1 the left single scroll, the centered single scroll and the left double scroll, plane S 2 is crossed by the right single attractor, the centered single scroll and the right double scroll. It is important to note that the double scroll atractors do not cross both planes, only the centered single-scroll crosses S 1 and S 2 , and none of the attractors share crossing points. This allows to guarantee the coexistence of each attractor.

Basin of Attaction Evolution due to the Fractional Order
To ilustrate the described behavior the corresponding basins of attraction are calculated, four of which are shown in Figure 7, each one considering different integration values. This allows us to see the deformation of the basin of attraction of the three-state system into a riddled basin of attraction [50], while as the order of integration increases, the double-scroll states become more probable until their regions of attraction collide and the global basin of the monostable three-state system emerges.
The basins of attraction shown in Figure 7 are based on the discretization of the attractors shown in Figure 4, in which all possible behaviors of the multi-scroll system are represented when modeled with fractional order derivatives. Each of the dynamics is assigned a number, 5 for the positive attractor with two scrolls, −5 for its mirror behavior. The dynamics with one scroll are assigned the values −3, 0, 3 depending on its position. The monostable attractor is associated with 10. In this way, and by assigning the colors of the attractors to the basin of attraction, the combinations of behaviors present in the system, which directly related to the order of integration used, are made visible in each basin of attraction.

Discussion
The analysis of the bifurcation diagrams presented in Figure 3, as well as the basins of attraction presented in Figure 7, show that the system can exhibit the coexistence of different attractors and that the mechanism leading from monostable to multistable dynamics is based on the reduction of the accessible domains in the trajectory of the system, leading to the fragmentation of the basin of attraction and, consequently, to the appearance of the different multistable states as presented in Figure 4.
From the recurrence points obtained from the Poincaré sections, shown in Figure 6, it can be seen that all the attractors shown in Figure 4 coexist, since they do not have the same recurrence points, which proves their uniqueness.
Due to the nature of the PWL function used in this work, it is possible to extend the number of domains and scrolls to n, and since the defined dynamics in each domain is identical, the existence of n + 2 coexisting attractors is guaranteed. To this end, consider a multi-scroll system with a 9-scroll attractor for its integer-order version, where the only difference from Eq. (10) is replacing a PWL function with more segmentations as described in Eq. (11) [51,52]. Under this assumption and using the same values as described in the previous sections, Figure 8 shows the multistable coexisting attractors and the global monostable dynamics for a 9-scroll system described by Eq. (9)(10)(11). The integer order attractor with 9 scrolls is shown in Figure 8(a), while the 9 coexisting single scroll dynamics are shown in Figure 8(b). The double scroll dynamics, coexisting in the left and right edges of the global dynamics are shown in Figure 8(c). Note that obtaining each attractor involves a fixed integration order q = 0.8963 and different initial conditions.

Conclusions
In the present paper, the effect of using fractional order derivatives in a PWL multi-scroll system was studied in detail. By analyzing the evolution of the attractor in the construction of bifurcation diagrams, it was possible to discover the route by which the system changes its global dynamics and moves from monostability to multistability. At the same time, it was discovered that such a path to multistability involves the presence of up to n + 2 attractors, rather than just n coexisting states as was previously known. It should be noted that n refers to the number of integer-order scrolls generated by the system. The coexistence of the n + 2 attractors is also demonstrated by the associated basins of attraction, which show the qualitative changes that the system undergoes due to fractional order.
The existence and uniqueness of the attractors is proved by means of the analysis of Poincaré sections, for each of the multistable attractors. These results confirm that the dynamics coexist in the different domains of the system, but they do not live in the same space, which is visualized by the location of the entry and exit points of the different domains of the system. It should be emphasized that the new coexisting states reported here have not been described in previous studies for both integer-order and fractional-order PWL systems.
To describe the mechanism of the occurrence of multistable behavior, one must assume that each of the equilibrium points can be considered as a potential well. In this sense, the reduction of the integration order leads to a reduction of the potential energy of the system, resulting in a segmentation of the basin of attraction into as many states as there are defined regions in the integer order system. At the same time, the occurrence of the double-scroll states at the ends of the attractor is due to the fact that there is no domain at the ends of the system. As a result, under certain initial conditions, shown in the basin of attraction, the dynamics can "gain"enough potential force to leave one domain and visit the next, but it does not have the necessary energy to restore the global dynamics of the system, resulting in a loop between two domains.
The corresponding analogical implementation to compare the presented results, as well as the extension of the results to another type of multiscroll systems, is proposed as future work. The results presented here extend the panorama of the effects of using fractional order derivatives and suggest the development of further and better studies to understand their effects on dynamical systems. realization of a post-doctoral stay at Dynamical Systems Laboratory, University of Guadalajara.

Declarations
The authors have no relevant financial or nonfinancial interests to disclose. All authors participated in the conception and design of the study. Material preparation, data collection, and analysis were performed by José Luis Echenausía Monroy. The first draft of the manuscript was written by JLEM, and all authors commented on earlier versions of the manuscript. All authors read and approved the final manuscript. Data sets generated and/or analyzed during the current study are available from the corresponding authors.