Explosive behaviors on coupled fractional-order system

Fractional derivatives provide a prominent platform for various chemical and physical system with memory and hereditary properties, while most of the previous differential systems used to describe dynamic phenomena including oscillation quenching are integer order. Here, effects of fractional derivative on the transition process from oscillatory state to stationary state are illustrated for the first time on mean-filed coupled oscillators. It is found the fractional derivative could induce the emergence of a first-order discrete transition with hysteresis between oscillatory and stationary state. However, if the fractional derivative is smaller than the critical value, the transition will be invertible. Besides, the theoretical conditions for the steady state are calculated via Lyapunov indirect method which probe that, the backward transition point is unrelated to mean-field density. Our result is a step forward in enlightening the control mechanism of explosive phenomenon, which is of great importance to highlight the function of fractional-order derivative in the emergence of collective behaviors on coupled nonlinear model.


Introduction
Spontaneous symmetry breaking in ensembles of interacting dynamic subsystems is an evergreen challenging problem on nonlinear science, which has widely applicability in diverse fields like physics, biology, and engineering [18,28,33,38]. In reality, spontaneous symmetry breaking issue in amplitude and phase variation on a coupled identical unit based on the breaking of the general translational and permutational symmetry of coupled oscillators, has induced several emergent collective phenomena such as chimera states [1,16,19], synchronization [12,30], extreme events [14,29] in natural and dynamic system. Recently, an ever-increasing interest on symmetry breaking phenomenon has been witnessed by the enormous burst studies on explosive synchronization (ES) and fractional differential equations.
The emergence of explosive synchronization (ES) on nonlinear dynamic system, which is characterized by the sudden (discontinuous) and irreversible transition between incoherence and coherence state as the variation of coupling strength, injects fresh blood to the research community to investigate the phase spontaneous symmetry breaking on coupled units. This kind of sharp and irreversible transition was first reported on phase model [13], where scale-free network and posi-tive correlation between the frequencies of the oscillators and their degrees were expressed explicitly as the two key factors of ES. Later, an important work [46] that generated a more general framework, frequencyweighted coupling where the before two ingredients are included as specific cases, was established. It has made ES become a topic of great interest among researchers on nonlinear science and fascinated plenty of significant studies [5,8,20,45] and experimental evidence [10,15,26,44]. Besides, a new phenomenon of explosive oscillation quenching on a frequency-weighted networked system [7], named as Explosive death (ED), grabs the attention of researchers owing to its rich background on complex systems [4,6,9]. In contrast to the classical continuous oscillation quenching transitions [32,48], ED implies the abrupt and irreversible process from oscillatory state to suppression of oscillation. Specifically, on account of the properties of the steady state, ED can be classified into explosive amplitude death and explosive oscillation death [17]. In fact, amplitude death (AD) appears when the coupled oscillators stabilize at a common homogeneous steady state (HSS), and when the stable HSS is associated with the coupling strength, nontrivial amplitude death (NAD) is observed. Whereas, oscillation death (OD) refers to that the system stays at the coupling related inhomogeneous steady state (IHSS), producing via the symmetry breaking, where all the oscillators occupy diverse stable branches of IHSS. Hence, when the system undergoes a first-order discontinuous transition from oscillatory to AD, explosive amplitude death occurs. While explosive oscillation death emerges when the stable state is OD. It is worth mentioning that, the rapid development of explosive behaviors has benefited a lot from the practical application on many complex networked systems and a great deal of research has been conducted for providing an effective control method on ED with different coupling schemes [11,23,34,[40][41][42]47].
As we all know, fractional differential equations play a vital role in many fields of science and engineering owing to the fact that the dynamics of many system can be rendered better by the fractional derivative [2,35,37]. In detail, most processes in real world with special materials and chemical properties are more likely to show fractional behaviors, such as properties of memory, non-locality and history-dependence. In the past few decades, an increasing number of efforts have been made to fractional-order system, especially in viscoelastic materials [3], fractional control prob-lems [22,31], relaxation processes in sensory adaptation [36] and so on. In fact, the extensive studies focusing on explosive phenomenon are done in integer coupling oscillators, while the effect of fractional derivative on the spontaneous symmetry breaking including ED has fundamental importance to the analysis and control of coupled system. So far, the contribution of fractional derivative on explosive phenomenon has not been revealed. Thus, in order to give a mechanistic insight to the explosive behavior, we elucidate a strategy of fractional-order mean-field coupling system in this paper that is outlines as follows. In Sect. 2, the fractional differential equations with Caputo derivative on mean-filed coupling system is elaborated detailly. In Sect. 3, the transition processes from oscillatory to stable steady state are presented numerically and theoretically. Furthermore, the roles of fractional derivative and mean-field density on explosive behaviors are revealed. Finally, the conclusion is summarized in Sect. 4.

Mathematical description
Without loss of generality, we analyze the van der Pol oscillation equation that was first proposed by Balthasar van der Pol [39] in vacuum tube circuits. In fact, the classical van der Pol oscillator described by the second-order differential equation:ẍ (t) − b 1 − x 2 (t) ẋ (t) + x (t) = 0 severs as a basic model of self-sustained oscillatory system, which is also known as unforced van der Pol equation. Recently, inspired by the development of fractional calculus, the fractional derivatives are applied in the van der Pol model in various forms. In particular, by introducing the capacitance via a fractance in the nonlinear RLC circuit, Pereira et al. [27] considered the fractional-order van der Pol (FVDP) as dt + x (t) = 0, 1 < α < 2 represents the fractional order. Meanwhile, another version of FVDP with a fractional damping term was obtained by Xie and Lin [43] as Subsequently, Mishra et al. [25] presented a general version of FVDP, namely, sequential fractional van der Pol equation, which is described as Further, the state space formulations are generated as, where D α t x i and D α t y i are the Caputo derivatives of x i and y i , 0 < α < 1 indicates α-order time derivative, respectively. In fact, the Caputo's definition of fractional derivative [21] is suitable to use in this paper on account of its superiority for the practice problem in engineering, which is defined as follows, where (·) is the Gamma function and n − 1 < α < n ∈ Z + . As the investigation of differential equations continued, the fractional-order counterpart of van der Pol oscillator has drawn more and more attention, many researchers in different branches have studied the dynamical properties of chaos and other complicated emergent collective phenomena.
Here, we intend to deal with the FVDP model in the form on Eq. (1). Owing to the fact that our main concern is whether explosive death could appear in the fractional-order dynamical system via a first-order discontinuous way, we consider a nonlinear model of N identical FVDP oscillators coupled mutually through a mean-field interaction. The dynamic is governed as follows, Herein, i = 1, 2, . . . N represents the number of coupled oscillators. The value Q is the intensity of the mean-filed with 0 ≤ Q < 1 and k (k > 0) indicates the coupling strength. The system parameter is b (b > 0) which shows the nonlinearity and damping strength of the single van der Pol oscillators and expresses a limit cycle.x = N i=1 x i /N purposed the mean-filed of the state variable x. Apparently, in the limit case α → 1 − , the first-order derivativeẋ can be obtained from Eq.
(2). That is,ẋ = lim α→1 − D α t x and Eq. (3) is reduced to the classical mean-field coupled van der Pol system, where oscillation quenching behaviors show with the variation of the coupling intensity.

Discontinuous transitions on coupled fractional-order system
We are interested here in the transition process from oscillatory state to the stationary state. To quantify the state at each coupling strength, i.e., to distinguish the oscillatory state and steady state, we compute a quantify indicator A (k) ∈ [0, 1] evaluated at a sample coupling strength k as the order parameter.
The variable a (k) indicates the difference between the global maximum and the minimum longtime average amplitude for a given coupling strength k when the system undergoes a sufficiently large transient time. The scaled average amplitude is given by, Thus, A (k) = 0 terms to the 'death' state and 0 < A (k) ≤ 1 represents the oscillatory state. It is noted that the order parameter A (k) is obtained numerically by the Runge-Kutta scheme. Owning to the fact that we are concentrating on investigating ED in fractionalorder system, we simulate the system adiabatically in both forward and backward directions through a continuation increase or decrease of coupling strength with step δk = 0.02. So far it has been reported [41] that in the limit case α → 1 − for N = 100 coupled van der Pol oscillators, the order parameter A (k) shows a sudden fall to zero with the increase of coupling strength. Still, as the coupling strength decrease, A (k) performs an abrupt jump from zero to a value which is repainted in Fig. 1a in the case of N = 2, indicating the firstorder transition appears in the integer order mean-filed coupling system is independent of the number of oscillators. On further observation of the time series near the critical transition points in Fig. 1b, this phenomenon is identified as explosive NAD. For the proposed consideration, we first consider the numerical results on a system of N = 2 cou- pled fractional-order system with a fractional derivative α = 0.8 in Fig. 2. As expected, we observe that the order parameter A (k) still exhibit a sharp, irreversible transition whenever we increase the coupling strength k from zero to the maximum values k max or decrease the k from k max to zero in Fig. 2a. Most certainly, these two sudden jumps emergent at different value of coupling strength k, showing a hysteresis area (HA) encircled by the forward and backward transition points k f , k b . However, when we reduce the coupling strength from the maximum k stacked in the forward direction, the transition path in the backward direction (blue lines) does not overlap the changing curve in the forward continuation (orange lines) after going through the backward transition point. That is, the oscillators cannot be restored to the initial state when the fractional derivative α = 0.8 was introduced, which is definitely a far cry from the previous research. Note that both the critical forward (the yellow star) and backward (the green cross) transition points are moving advanced compared to the integer order system in Fig. 1a. To better characterize the properties of the emerging explosive death, we then portray the corresponding time transient of the coupled system near the critical forward and backward transition points in Fig. 2b1 and b2. We can infer that the oscillators synchronize for a smaller value of k. Then with the coupling strength getting larger, the synchronous behavior disappears immediately after the critical forward transition point k f = 1.0. Interestingly, at the same time, the oscillators stabilized at two inhomogeneous steady branches. Namely, once the fractional derivative exists, the system will preferentially stabilize at the OD state rather than the homogenous equilibrium. In fact, when the coupling strength k is becoming smaller with the decrement, Fig. 2b2 reveals that the OD state is lost as soon as the backward transition point k b = 0.72 (less than k f ) reaches. It can be clearly seen that the system oscillates with two opposite amplitudes after this critical coupling strength, indicating that the system reaches an anti-synchrony state rather than initial synchronization. This is not surprising, since we have determined that the oscillators are stable on two branches in the forward transition, we further choose the inhomogeneous steady state as the initial state of k = k max in the backward direction rather than random selection like forward process, thus the anti-synchronized state is easier to obtain. Thus, not only the hysteresis area of OD and synchronized state is observed, but also a fresh region of synchronization and anti-synchronization appears. Moreover, the results above are independent of the number of coupled oscillators. In fact, in the case of N = 100, the explosive OD appears in both system and the forward transition points are the same as the case N = 2. Meanwhile, in the case of N = 2 , the two oscillators stable at two branches that are symmetric about the origin. However, oscillators distribute in the upper and lower branches which is asymmetric about the origin for N = 100 . Hence, the initial state of N = 2 and N = 100 in the backward direction is different, inducing the diverse backward transition points. Besides, the system experiences a sudden and irreversible transition from oscillatory state to chimera death state, which has a great difference from the case for N = 2 (see details in Sect. 3.3).

Amplification of the hysteresis area
When the explosive discontinuous transition happens, there must exist divergent critical transition points in two directions owing to the adiabatic methods. The region (nondimensional area) enclosed by these two points is so-called hysteresis area which is known as the coexistence of OD and synchronization in the previous section. To achieve more enhancement with the characteristics of HA in the fractional-order system, we perform the variation of forward transition point k f , backward transition point k b and HA width along with the decrease of fractional derivative α from 1 to 0 by decrement δα = 0.01 in Fig. 3 for a fixed Q = 0.5, and some interesting conclusions are drawn. Specifically, a tall hysteresis area is clearly found to pop up for α = 0.99 in Fig. 3a, meaning that the explosive OD happens even for α = 0.99. In other words, there is bound to exhibit explosive OD not explosive NAD once the fractional derivative is introduced. Then we decrease α and observe a uniform reduction of the hysteresis region. Until the threshold value α min = 0.72 arrive, the forward and backward transition points coincide, indicating that the appearance of hysteresis loop made up of OD and oscillatory state is closely connected with fractional derivatives. In fact, the critical forward and backward transition points in the current paper are described as the points after the state variation, thus, when k b − k f = δk = 0.02 is established, we believe the two critical points overlap. It is tentatively suggested that, the hysteresis region is decreasing with the reduce of fractional derivative. More to the point, the fractional derivative α > 0.72 is what we see as a qualifying condition to ensure the occurrence of explosive OD. Besides, we note that, with the increase of α, the hysteresis area width enlarges linearly. Thereby, the variation of HA width with the fractional derivative is well fitted by the linear function H A width = 3.2699α − 2.3277 (dotted line in Fig. 3b).
Encouraged by the interesting observation on the occurred condition of the explosive OD in fractionalorder system, we now proceed to vary the intensity of the mean-filed Q and investigate the peculiarity of the explosive transition. In fact, we start with the value Q = 0 and increase it with increment δ Q = 0.1 and go up to Q = 0.9 as described in Fig. 4. We note that, the threshold value α min doesn't growth with the enlargement of the mean-field intensity Q. Further, the intensity of the mean-filed could implement as a regulator to achieve the reviving of explosive behavior. In other words, to take a small example, in the case Q = 0.1, the explosive OD will vanish if the fractional derivative fit with α < 0.8. Then we increase Q to Q = 0.4, the fractional derivative α in the interval of 0.72 < α ≤ 0.8 could induce explosive behavior. In this way, the value Q could adjust the occurrence of the first-order transition. Thus, it is verified that the mean-filed intensity plays a vital role of controlling the qualifying condition of the explosive OD.
Moreover, for detail investigation, we then compute the critical transition points during the forward and backward continuation with the increase of the meanfiled intensity for different fractional derivative α in In fact, it is not surprising to note that the HA width in Fig. 5c performs the similar results with the forward transition points since it is obtained by subtracting the value of backward transition point from the value of forward transition point. Meaning that the fractional derivative and the mean-field intensity could influence the appearance of the discontinuous transition from synchronization to stationary state in the forward direction. Besides, the fractional derivative also acts an important role to adjust the backward transition from OD to anti-synchronized state. However, the backward transition point is entirely unrelated to the mean-filed intensity that is in good agreement with the theoretical results (see details in Sect. 3.2), which indicates that the backward transition point is independent of the mean-filed intensity Q.

Analytical analysis of critical transition point
To further enunciate the transition in the explosive behavior, we move to some theoretical analysis via Lyapunov indirect method in the next step of our study. Specifically, the fixed points of the two coupled fractional oscillators can be divided into three categories: (a) trivial steady state E 0 = (0, 0, 0, 0). x * = 1 − 1 bk and y * = kx * , which is appearing due to the symmetry breaking. Note that the IHSS exist under the condition k > 1/b . Moreover, the IHSS is independent of the mean-field density, which is showing detailly from the form of the inhomogeneous fixed point. From the results presented above, it is clear observed that the stationary state of meanfield coupling fractional-order system is OD. Thus, we mainly adopt the Lyapunov indirect method on the coupling related IHSS E * . Thereby, we perform the linear stability analysis in the inhomogeneous equilibrium E * to obtain the stable boundary, in other words, the critical transition point. The characteristic equation is given by, Then the corresponding eigenvalues are described as, Herein, According to the stability theorem of fractional systems proposed by Matignon [24], the fixed points E * is asymptotically stable if all the eigenvalues λ * m (m = 1, 2, 3, 4) corresponding to the linear matrix content the criterion, That is, if the absolute values of all the eigenvalues corresponding to the linear matrix are greater than απ/2 , the nonlinear system is asymptotically stable. Hence, the stable inhomogeneous equilibrium of coupled fractional-order oscillators, i.e., OD state can be determined as following, . Therefore, an explicit illustration of the backward critical transition point is exhibited. In fact, when the damping coefficient b is smaller than √ 1 − Q, the backward transition point k b is 1/b . Otherwise if the damping coefficient b and the mean-field density Q meet with the situation b > √ 1 − Q, a relationship between the fractional derivative α and the critical backward transition point k b is drawing clearly by the expression in Eq. (8). That is, the critical points in the backward direction should match the condition α < One final thing to note here is that, on the condition of Figs. 3 and 5, we can examine that if the backward critical transition point satisfies the crite- that is marked as black solid line in Fig. 3a, the oscillators are stable. In addition, we observe that the stable condition is independent of the mean-flied density. However, the backward transition points can be divided into two types via the relationship of b and √ 1 − Q. Besides, one can find that the criteria have no relationship to the mean-field density Q which is fitted with the numerical results in Fig. 5b. Therefore, it verified that the backward transition point is independent of the mean-field intensity in the current paper.

Emergence of explosive chimera death
Based on the above observation, explosive NAD exist in integer order system really does not appear in the fractional mean-flied coupling oscillators, but oppositely a sudden and irreversible transition from oscillatory to OD state once the fractional derivative is introduced. However, the system reported in the previous sections is only a special case for N = 2, it is valuable to study how the system changes with the degrees-of-freedom of the coupled oscillators. For simplicity, we divided the network size into odd and even cases and consider the transition processes in Fig. 6, where we take N = 100 and N = 95 as examples. Through observing the variation of order parameter and time transient near the transition points, we find that, the sudden and irreversible transition from synchronous to OD state happens in both odd and even cases at the forward process. And in the backward direction, with the decrease of the coupling strength the system undergoes a sharp jump from OD to synchronization. Meanwhile, a hysteresis area of synchronization and OD coexistence. Besides, the critical forward transition points are the same in both cases, while the backward critical transition points have a fine distinction between N = 100 and N = 95. It is interesting to note that, in the backward direction, oscillators go back to synchronous oscillation not anti-synchronization showing in Fig. 2 in the backward direction, which has a great difference from the transition for N = 2. Therefore, the coexistence state of synchronization and anti-synchronization appearing in Sect. 3.1 will not exist for high dimensional system.
To provide more insight into the discontinuous firstorder transition appearing in the high dimensional fractional system, we examine the space-time plots and the snapshots of the variable X i near the forward and backward transition point in Figs. 7 and 8. Still the in-phase synchronized state is well pronounced in the oscillatory state for both odd and even cases. Interestingly, through further analyzing the space-time and snapshots at the stable inhomogeneous steady state, we observe that the N identical oscillators split into two subpopulations: the neighboring oscillators distribute Fig. 7 a-b Space-time plots for the variable X i (t) near the forward critical transition points for different network size a1-b1 N = 100 and a2-b2 N = 95. c Snapshots at t = 1000 in the forward direction to show the explosive chimera death the same branch of inhomogeneous steady state on one domain, while the neighboring elements are populating randomly on the two branches of the inhomogeneous steady state in the other subpopulations. Hence, the stationary state is characterized by the coexistence of spatially coherent oscillation death and spatially incoherent oscillation death, which is the typical characteristic of chimera states. In other words, the stationary state is further identified as chimera death (CD) via Figs. 7 and 8. Therefore, the oscillators undergo the sudden jump between synchronized state and CD, while a hysteresis area of CD and synchronization exist, that is, the explosive chimera death emerges in N coupled fractional oscillators. In more detail, for high dimensional cases N = 100 and N = 95, the fractional system experiences a first-order transition from in-phase synchronized state to chimera death state, which has a great difference from the case for N = 2. Moreover, although the explosive CD exist in both even and odd situations, the distribution of nodes which is stable at coherent and incoherent oscillation death is completely random. When we reduce the coupling strength from k max to zero, we choose the final state of the forward transition rather than the uniformly symmetrical state in Sect. 3.1 as the initial state. Thus, there are subtle differences in the backward critical transition points for N = 100 and N = 95.
For detail investigation of the explosive transition evolution with the number of network node N , we examine the variation of forward, backward transition points as well as the hysteresis area width with the increase of network size in Fig. 9. Interestingly, the cru-  cial transition points are independent of the parity of the number of oscillators to some extent. In particular, the forward transition points even remain unchanged for any network size we showed. Whereas, as for the stationary state in the forward direction, the distribution of the coherent OD and incoherent OD for diverse network size is random, which is the initial state of the backward process. Thus, the backward transition points fluctuate slightly in a range. Meanwhile, the width of the hysteresis area composed of synchronized state and chimera death state also performs floating in a small interval as showing in Fig. 9b.

Conclusion
The present study has dealt with the problem of explosive phenomenon on a system of mean-filed coupling fractional-order nonlinear oscillators. We demonstrate the setting of a sudden and irreversible transition from oscillatory to stationary states, where it is observable even when the fractional derivative is proposed to describe the chemical properties and special materials in real world. It is uncovered that, once the fractional derivative is introduced, the system prefers to stabilize at the inhomogeneous steady state rather than the homogenous equilibrium. Meanwhile, it is indicated that fractional derivative exhibits an amplification of the hysteresis area enclosed by the critical forward and backward transition points of the explosive transition. Indeed, our result is not limit to the case of two coupled oscillators, which is fully robust against the large amount of the oscillators. In more detail, for high dimensional system, oscillators undergo the abrupt transition from synchronized state to chimera death state and result in a coexistence state of synchronization and CD state, which is independent of the parity of the number of oscillators. Moreover, for further consideration on the role of the fractional derivative, it could be revealed that, there is a minimum fractional derivative for inducing explosive behavior when the mean-filed density is fixed. In other words, the appearance of the first-order discontinuous transition could be adjusted via the fractional derivative. Interestingly, the threshold value α min does not increase with the enhancement of the mean-field density Q. Thus, an appropriate mean-filed density and sufficiently large fractional derivative is liable to induce the emergence of explosive behavior.
Besides, it is further verified that the backward transition point has no relationship with the mean-filed density numerically and theoretically. The mean-filed density could just regulate the forward direction to control the coexistence area of oscillatory and stationary state. Notably, our findings mainly focus on a new attempt in the fields of investigation and application of the explosive transition in natural systems showing fractional behaviors. Our results about understanding the role of fractional derivative on collective behaviors especially in oscillation quenching can be constructive in absorbing the underlying mechanism of many realistic chemical and physical systems with properties of memory, non-locality and history dependence. tation of Northwestern Polytechnical University (Grant Nos. CX2021035). The authors would like to thank the anonymous referees for their efforts and valuable comments.
Funding National Natural Science Foundation of China (Grant Nos. 11772254, 11972288) and Innovation Foundation for Doctor Dissertation of Northwestern Polytechnical University (Grant Nos. CX2021035) Data availability Data will be made available on reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest, with respected to the research, authorship and publication.