Non-Smooth Dynamic Behaviors as well as the Generation Mechanisms in a Modiﬁed Filippov-Type Chua’s Circuit with a Low-Frequency External Excitation

: The main purpose of this paper is to study point-cycle type bistability as well as induced periodic bursting oscillations by taking a modiﬁed Filippov-type Chua’s circuit system with a low-frequency external excitation as an example. Two different kinds of bistable structures in the fast subsystem are obtained via conventional bifurcation analyses; meanwhile, nonconventional bifurcations are also employed to explain the nonsmooth structures in the bistability. In the following numerical investigations, dynamic evolutions of the full system are presented by regarding the excitation amplitude and frequency as analysis parameters. As a consequence, we can ﬁnd that the classiﬁcation method for periodic bursting oscillations in smooth systems is not completely applicable when nonconventional bifurcations such as the sliding bifurcations and persistence bifurcation are involved; in addition, it should be pointed out that the emergence of the bursting oscillation does not completely depend on bifurcations under the point-cycle bistable structure in this paper. It is predicted that there may be other unrevealed slow–fast transition mechanisms worthy of further study.


Introduction
Many practical problems, such as neuron activities [1][2][3], surface oxidation reactions and autocatalytic chemical reactions [4][5][6], memristor-based circuits [7,8], and so on, can often be described as slow-fast dynamical systems (SFDSs), written in the generalized form where · = d dτ , u ∈ R N are the fast variables, v ∈ R M are the slow variables, the positive parameter 0 < ε 1 measures the timescale difference between u and v, and f : R N × R M × R → R N and g : R N × R M × R → R N are the smooth functions.
Compared with single-timescale dynamical systems, one can find that the system responses of (1) may be much more complicated, such as in the bursting oscillation patterns concerned in this paper. A bursting oscillation generally refers to the oscillation mode in which the trajectories of the whole system alternate between the quiescent states(QSs) and the spiking states (SPs) periodically. QS means that the trajectory converges to a stable equilibrium, while SP means that the trajectory oscillates along with a limit cycle. Since such complicated oscillation modes were reported in a large number of experiments and numerical simulations, the related dynamical evolutions as well as the generation mechanisms confused the researchers for a long time based on the fact that the mature analysis techniques in single-timescale dynamical systems cannot be directly applied to (1) when bursting phenomena are present.
Fortunately, an incredibly useful approach based on timescale separation, i.e., the classical slow-fast analysis method, was introduced by Rinzel to expound the dynamics underlying bursting in neurons and pancreatic β-cells [9][10][11]. It is feasible to understand a burster by employing nonlinear dynamics techniques and bifurcation theory. Taking system (1) as an example, potential bursting oscillations may be well explained by executing the following two stages: firstly, dividing the full system (1) into a fast subsysteṁ u = f(u, v, ε) and a slow subsystemv = εg(u, v, ε) via timescale separation; secondly, analyzing the attractors as well as the bifurcations in the fast subsystem by regarding the slow variables v as bifurcation parameters. Then, it can be found that a potential burster in System (1) may be well interpreted as the dynamical evolution process in which the trajectory modulated by a slow subsystem periodically alternates between a spiking attractor and the following quiescent attractor (for instance, from a stable limit cycle to a stable equilibrium regime, and so forth), which is caused by bifurcations in the fast subsystem. By adopting such a method, some extensive generation mechanisms of bursters have been listed by Izhikevich in his works [12,13]. More importantly, he also proposes that a burster can be named by using the two important bifurcations leading to the transitions between SPs and QSs, which can effectively distinguish and classify different bursters and their generation mechanisms and may provide great convenience for studies related to bursting phenomena.
Along with the successional exploration of various bursting forms and the generation mechanisms in slow-fast dynamics, accumulating evidence suggests that bistability, referring to when one attractor coexists with another attractor under certain parameter conditions, may play an important role in the dynamical evolutions of bursters [14,15]. The famous binocular rivalry, which reflects the human event for which a person reports an alternation between two competing percepts as opposed to a mixture of them when his/her eyes are exposed simultaneously to two significantly different images, can be explained theoretically in the view of dynamics by the so-called winner-take-all bursting transition mechanism, in which the trajectory may converge to one of the bistable equilibrium branches to behave in QS with different parameter conditions [16,17]. A novel bursting form named mixed bursting oscillation has been observed in [18], which is characterized by the fact that different SP patterns can be performed in a period bursting oscillation. By using slow-fast decomposition, Duan et al. pointed out that different bistability structures may lead to different mixed bursting patterns [19].
Note that most of the above results are verified in continuous slow-fast dynamical systems; however, some systematic mutations may also appear in slow-fast dynamical systems, for instance, nonsmooth stillness in nonlinear energy sinks [20], dry friction in an oscillator with low-frequency excitation [21], and switch-like interactions in gene regulatory networks [22]. Then, corresponding mathematical models are discontinuous, and the conventional bifurcation theory is nearly helpless for the understanding of the slow-fast dynamics as well as the generation mechanism.
Particularly in today's third generation of neural networks (referring to spiking neural networks), although the adoption of piecewise smooth integrate-and-fire models for a single neuron bring great convenience to the computability and application of neural networks (for more details, see [23] and the references therein), the reset motion in neuron models and gap-junction coupling among neurons perform intricately discontinuous and implicit nonlinear mechanisms.
Obviously, a deep understanding of bursting phenomena as well as the mechanism in discontinuous slow-fast dynamical systems is very meaningful to slow-fast dynamics as well as their practical applications. Based on this purpose, we focus on the slow-fast dynamical systems with discontinuous vector fields (which can also be called Filippov-type slow-fast dynamical systems). The standard mathematical model can be expressed bẏ in which, the fast subsystem possesses a discontinuous boundary (switching boundary) H(u) = 0 that divides the fast subsystem into two smooth subsystems. Apparently, when contacting with H(u) = 0, attractors in the fast subsystem of (2) may undergo nonconventional bifurcations, leading to System (2) exhibiting unique nonsmooth dynamic behaviors. This paper aims to discuss those dynamic behaviors as well as the generation mechanisms in (2). The remainder of this paper is organized as follows. In Section 2, a Filippov-type slow-fast dynamical system is presented via establishing a modified Chua's circuit with a low-frequency external excitation. Necessary bifurcation analyses are presented in Section 3. Based on Section 3, the attractors of the fast subsystem and their stabilities are discussed with property parameter settings. Various dynamical behaviors observed in the full system as well as the underlying generation mechanisms are discussed via numerical investigations in Section 4. We draw some conclusions and discussions in Section 5.

Mathematical Model
The classical Chua's circuit, designed by Chua, composed of two capacitors, an inductor and a Chua's diode, is a famous 3D nonlinear electronic circuit to show chaotic behavior [24,25]. The circuit as well as the various modifications are often taken as examples in the research on slow-fast dynamics. For instance, canards and chaotic bursting can be revealed in memristor-based Chua's circuits [26], or more to the point, those modified Chua's circuits with low-frequency excitations are employed for the investigation of various bursting oscillations as well as the generation mechanisms via slow-fast decomposition [27][28][29].
To reveal the influence of non-smoothness on the dynamics with two scales, here we consider one modified Chua's circuit system by introducing a low-frequency current source and a nonlinear resistor R N with piecewise characteristics, shown in Figure 1. The related mathematical mode can be written as where i S = I m sin(ω t) denotes a sinusoidal AC power source, which represents an external low-frequency excitation, and R N is the nonlinear resistor designed in Figure 2. Figure 2. The design of the nonlinear resistor R N in Figure 1.
Let i be the current flowing through the resistor R N ; the v − i characteristic can be written as By employing the transformations τ = 1 Assuming that W sin(Ω t) is considered a low-frequency excitation, i.e., 0 < Ω 1, indicating that the whole excitation changes periodically on a slow timescale in [−W, +W], System (5), which is one of the slow-fast dynamical systems with two timescales coupling in the frequency domain, can also be called a Filippov-type slow-fast dynamical system.
According to the generalized slow-fast analysis method [29], the attractors and their bifurcations parameterized by the slow variable w = W sin(Ω t) are the core to understand potential bursting oscillations in (5). However, the discontinuity boundary H(X) = 0 (labeled as Σ) divides the state space of the fast subsystem, denoted as Π(X), into two smooth sub-regions: R − := {X|x < 0} and R + := {X|x > 0}, i.e., Π = R − ∪ R + ∪ Σ. At the same time, in the smooth sub-regions R − and R + , there are two smooth sub-systems S − and S + , expressed by and The corresponding vector fields are denoted as F − and F + , respectively. Meanwhile, for a point X Σ ∈ Σ, using the directional derivatives of H(X) with respect to the two vector fields on both sides of Σ, expressed as ∂X , where ·, · is the vector inner product), we have the following results on Σ with the assumption that δ < 0.
When H, F + H, F − > 0, there exist two sewing regions where the trajectory may behave in a crossing motion after it contacts with Σ + or Σ − ; when H, F + H, F − < 0, there exists a sliding region Moreover, by introducing an auxiliary parameter Q ∈ [−1, 1] via Utkin's equivalent control method [30], one 2D sliding vector field in Σ S can be computed as corresponding to the sliding subsystem Particularly, Q = −1 and Q = 1, respectively, correspond to two important sliding boundaries For one bursting oscillation in System (5), when the trajectory is in the smooth subregion R − (or R + ), it will be driven by the smooth subsystem S − (or S + ), and its dynamics evolution process can be explained by the smooth dynamics theory. However, when the trajectory interacts with the non-smooth interface Σ, unique nonsmooth dynamics in Filippov systems, including sliding motions as well as nonconventional bifurcations, may appear and influence the bursting structure.

Stability Analysis of the Slow Manifolds
As the vector fields of the fast subsystems (6) and (7) are origin-symmetry, the equilibrium point can be labelled in a unified form X 0 (x 0 , 0, −x 0 ). Then, the slow manifold is written as which consists of three parts: 1.
The equilibrium points of S − and S + in the smooth region are called admissible equilibrium points, denoted as AE − and AE + , respectively, which can be expressed as 2.
The equilibrium of the sliding subsystem (11) located in the sliding region is called a pseudo-equilibrium point, denoted by PE, which can be expressed as 3.
The equilibrium points on the sliding boundary ∂Σ S − and ∂Σ S + are called boundary equilibria, denoted by BE − and BE + , respectively, which can been expressed as Note that the PE in the sliding region is equivalent to the equilibrium of the sliding subsystem, so the stabilities of AE ± and PE can be characterized by the corresponding characteristic equations.
For PE, the corresponding characteristic equation can be written as The eigenvalues of (18) are easy to derive as λ BE 1,2 = −0.5 ± 0.5 1 − 4 η. It is not difficult to find that the pseudo-equilibrium is always stable when η > 0.

Hopf Bifurcations on AE
For a stable AE on a slow manifold, labelled as AE := (x 0A , 0, −x 0A ; w), one important route to lose its stability is the occurrence of a Hopf bifurcation characterized by a pair of conjugate complex roots of (17) λ AE 1,2 = ρ(w) ± β(w) I crossing through the imaginary axis and the rest root λ AE 3 < 0. By omitting redundant calculations, the exact conditions for Hopf bifurcations on AE are given in the following.
The first equation guarantees (17) has a pair of pure imaginary roots, the second equation indicates the transversality, and the last one represents the nondegeneracy. Furthermore, if L 0 < 0, it is a supercritical Hopf bifurcation; if L 0 > 0, it is a subcritical Hopf bifurcation.

Nonconventional Bifurcations of the Slow Manifold
The two boundary equilibria BE ± (X 0B ; w B ) = (0, 0, 0; ∓δ) may connect PE with AE when w = ∓δ. Then, nonconventional bifurcations on the slow manifold, also called boundary equilibrium point bifurcations (BEBs), may be observed. Now, we take BE − as an example to exhibit more details.
At the boundary equilibrium BE − , the condition for the codimension-1 bifurcations of the boundary equilibrium must be satisfied in which the first is a non-degenerate condition to ensure that BE − is isolated in the vector field, while the second is a transversal condition. In this case, BE − connects AE − and PE, indicating that there is an admissible equilibrium ( and a pseudo-equilibrium point (X P− ; w * ) ∈ PE, which can be represented by introducing an auxiliary parameter q = (1−Q) Linearizing (21) and (22) at BE − (X * ; w * ), we have Since 2 δ < 0, the admissible equilibrium and the pseudo-equilibrium may collide into each other at BE − and disappear together, indicating the appearance of a non-smooth fold bifurcation; if α (1+a) 2 δ > 0, the admissible equilibrium and the pseudoequilibrium are located on the opposite side of BE − and transform to each other through BE − , which means the appearance of a persistence bifurcation. Thus far, Hopf bifurcations as well as BEBs on the slow manifold have been analyzed.

Identification and Localization of the Sliding Bifurcations
In the bursting oscillation, the spiking state usually corresponds to the process of the trajectory following the stable limit cycle, so the potential stable limit cycle can be regarded as the attractor of the spiking state. When the limit cycle interacts with the discontinuous interface, according to (10), the sliding vector field can be expressed as When H, F S = 0, we have Then, for the sliding region Σ S and sliding boundary ∂Σ S± on the discontinuity boundary, the following conclusions are obtained If we extend the vector field on the interface as It is obvious that Q ∈ [−1, 1] corresponds to the above conclusion, and further, in the nonsmooth crossing region, we have without losing generality, taking one stable limit cycle interacting with Σ as an example, there may exist two important contacting points: one is the point at which the trajectory comes into contact with Σ from the smooth region R + , while the other is the point at which the trajectory comes into contact with Σ from the smooth region R − , respectively corresponding to two Q values, denoted as Q re+ and Q re− . Obviously, one may find that the nonsmooth oscillation structures of the corresponding stable limit cycle may be directly reflected by the values of Q re+ and Q re− , and moreover, two curves formed by Q re+ and Q re− in the (w, Q) plane may also draw the nonsmooth dynamics evolution process via combining Q = ±1, indicating that sliding bifurcations can also be located and identified by Q re+ and Q re− . Such a numerical method can be called a nonsmooth returning map since it represents the first interacting points of trajectories of a limit cycle and discontinuity boundary by discrete points in (w, Q), and more details about this will be presented in numerical simulations.

Numerical Simulations
In this part, we turn to numerical study by using the method of the fourth-order Runge-Kutta. According to the theory analyses, when setting α = 1, a = 0.2, b = −1.6, c = 1, δ = −0.1, η = 0.1, obvious bistability can be observed in the fast subsystem, as shown in Figure 3.
It can be found that the four subcritical Hopf bifurcations, computed as subH 2± (w, x) = (±0.3627, ±0.2357) and subH 1± (w, x) = (±0.6429, ±0.9510), divide the slow manifold M into three parts: stable AE 1± , stable AE 3± , and unstable AE 2± . Meanwhile, two persistence bifurcations PB ± , located at BE ± (w, x) = (±0.1, 0), connect the stable pseudo-equilibrium branch AE 3± , respectively. Note that one spiking attractor, i.e., a stable limit cycle LC S represented by two extremes LC max S and LC min S , exists between w ≈ −0.6552 and w ≈ 0.6552. Then, the bistability in Figure 3a will be separated by the unstable limit cycles bifurcated from the subcritical Hopf bifurcations subH i± (i = 1, 2). Accordingly, the unstable limit cycles can be approximately found by the dichotomy between an extreme point of the stable limit cycle and the stable equilibrium point based on attractor basin theory, respectively labelled as LC U1± and LC U . Furthermore, the two fold bifurcations of limit cycles LPC ± occur via collisions between the spiking attractor LC S and the unstable limit cycles LC U1± bifurcated from HB 1± .
Meanwhile, based on the fact that the limit cycle LC S always contacts with Σ, we draw two Q curves Q re+ (light gray dashed line) and Q re− (red dashed line) in Figure 3b. The two curves Q re± interact with sliding boundaries ∂Σ S± at (w, Q) ≈ (±0.6161, ∓1), (w, Q) ≈ (±0.5307, ∓1) and (w, Q) ≈ (±0.2915, ∓1), and obvious variations of nonsmooth structures may be observed via the Q re± values when LC S passes through these points along with changing w, indicating that unique sliding bifurcations may appear. We now take the LC S located in w ∈ [0, 0.6552) as an example to give more details based on the geometry structures of sliding bifurcations in [30].
When w ∈ [0, 0.2951), both Q re+ and Q re− are located in the regions with |Q| > 1, i.e. sewing regions, then LC S may directly cross through Σ after it interacts with Σ, behaving in double crossing nonsmooth oscillation mode in a complete oscillation period, as illustrated by LC S with w = 0 in Figure 4a. Along with w increases over w = 0.2951, Q re− may enter into sliding region represented by |Q| < 1 via crossing through ∂Σ S+ (Q = 1) while Q re+ is still in the region Σ − (Q < −1), indicating that one sliding motion amd one crossing motion can be observed in LC S , as shown in Figure 4b with w = 0.38. Obviously, the nonsmooth oscillation mode has transformed to sliding crossing oscillation mode from double crossing oscillation mode, meaning that one sliding bifurcation has occurred on sliding boundary ∂Σ S+ when w = 0.2951. In such a sliding bifurcation structure, the intersection of Q re− and ∂Σ S+ actually represents the geometry structures that intersecting points of trajectories starting in one smooth subsystem and discontinuity boundary may continuously pass through the sliding boundary belonging to the opposite smooth subsystem, referring to the nonconventional bifurcation called crossing-sliding bifurcation. Therefore, we can say that LC S goes through crossing-sliding bifurcation when w = 0.2951. Similar bifurcation structure can also be observed at w = 0.5307, i.e., crossingsliding bifurcation of LC S also takes place when w continuously increases over w = 0.5307, leading to that LC S may behave in double sliding oscillation mode from crossing sliding oscillation mode after w enters into w ∈ [0.5307, 0.6161) since both Q re+ and Q re− are located in sliding region, as shown in Figure 4c with w = 0.58.
Particularly after w increases over w = 0.6161, Q re− may disappear on sliding boundary ∂Σ S− while Q re+ is still in the sliding region. Then, two sliding motions merge into one via the geometry structures for which intersecting points of trajectories starting in one smooth subsystem and discontinuity boundary may continuously approach the sliding boundary belonging to the same smooth subsystem and disappear ultimately, thus actually corresponding to the geometry structures of multi-sliding bifurcation [31]. Therefore, we may say that LC S goes through multi-sliding bifurcation when w = 0.6161, leading to the transition from double sliding oscillation mode to single sliding oscillation mode, as shown in Figure 4d with w = 0.63.
According to the symmetry, sliding bifurcations of LC S , which are located by the six intersecting points of Q re± and sliding boundaries, have been distinguished via the correspondence between the local structures of Q re± interacting with sliding boundaries and geometry structures of sliding bifurcations, including two multi-sliding bifurcations MS ± with w = ±0.6161 as well as four crossing sliding bifurcations CS 1± with w = ±0.5307 and CS 2± with w = ±0.2951, as shown in Figure 3b.
After clarifying the attractors and bifurcations of the fast subsystem under the given parameters, it can be found that there are three stable parts on the slow manifold M, which are further expressed as follows: Then, the stable slow manifold M iS (i = 1, 2, 3) and the spiking attractor LC S are divided into two monostability zones and three bistability zones by four sub-Hopf bifurcations HB i± (i = 1, 2) and two fold bifurcations of limit cycles LPC ± ; more details are given in Table 1.

Influence of Excitation Amplitude
In order to study the responses of the Filippov-type slow-fast dynamical system (5) with changing excitation amplitude W, we firstly set Ω = 0.001 1 and W ∈ (0, 1). According to the bifurcations in Figure 3 as well as the stability structures in Table 1, it can be found that there are four different parameter regions: W ∈ (0.6552, 1), W ∈ (0.3627, 0.6552], W ∈ (0.1, 0.3627] and W ∈ (0, 0.1]. Without loss of generality, we take the following four visiting modes of slow variable mode A: W = 0.8, mode B: W = 0.655, mode C: W = 0.3 and mode D: W = 0.1 as examples to exhibit more details. The corresponding responses of the waveforms of the full system are given in Figure 5. When W = 0.8, as shown in Figure 5a, one typical periodic bursting oscillation represented by the alternations between two quiescent states QS i (i = 1, 2) and two spiking states SP i (i = 1, 2) can be observed in the waveform. Based on the generalized slow-fast analysis method, Figure 6 gives the overlapping of the bursting oscillation and the one parameter bifurcation diagram in the (w = 0.8 sin(Ωτ), x) plane. It is not difficult to find that the two spiking states SP i (i = 1, 2) and the two quiescent states QS i (i = 1, 2) are formed by the trajectory, respectively, moving along with the spiking attractor LC S and stable admissible equilibrium branches M iS (i = 1, 3) on the slow manifold.
Furthermore, Figure 6b gives the transition mechanism from SP 2 to QS 1 and then to SP 1 , where the fold bifurcation of limit cycles LPC + leads to the vanishing of the spiking attractor LC S with w ≈ 0.6552 via colliding with the unstable limit cycle LC U1+ , further resulting in the performance of the transition from SP 2 to QS 1 . While passing through the subcritical Hopf bifurcation point HB 1+ with w ≈ 0.6429, the trajectory may gradually tend to the outside spiking attractor LC S after a common slow passage effect [25], since the slow manifold looses its stability via HB 1+ , resulting in the transition from QS 1 to SP 1 . According to that, one periodic bursting oscillation can be named and classified by the bifurcation that leads to the transition from quiescent state to the following spiking state as well as the bifurcation that leads to the end of the spiking state [12]. Here, the periodic bursting oscillation in Figure 5a can be called a symmetrical sub-Hopf/LPC burster. Since the spiking states SP i (i = 1, 2) are formed by the trajectories oscillating along with the spiking attractor LC S , sliding bifurcations presented in Figure 3b unavoidably affect the oscillation structures of the two spiking states, leading to the fact that the four nonsmooth oscillation modes in Figure 4 may be observed, as shown in Figure 6c. By drawing the returning maps Q re± of spiking state SP 1 , the spiking state SP 1 actually behaves in unique mixed oscillation modes via performing nonsmooth oscillation mode series: single sliding −→ double sliding −→ crossing sliding −→ double sliding −→ crossing sliding −→ double sliding −→ single sliding.
When W ≤ 0.6552, the two fold bifurcations of limit cycles LPC ± are no longer in the coverage area of w, resulting in the fact that the transition from LC S to M iS (i = 1, 3) cannot appear. As shown in Figure 5b, Figure 5c1 with initial values X in = (2, 2, 2) and Figure 5d1 with initial values X in = (2, 2, 2), only oscillations restricted on the stable cycle manifold can be observed in the last three cases once the full system trajectory converges to the spiking attractor LC S . Those oscillations are actually quasi-periodic oscillations that can be well proved by the closed Poincare maps via taking τ = 2 k π + π/2 (k ∈ N * ) as the section (see Figure 7).
Particularly, when W ≤ 0.3627, only the bistability consisting of M 2S and LC S lays in the convergence area of w, indicating that coexistence phenomena may be observed in the full system. For instance, as seen in Figure 5c1,c2 with W = 0.3, when setting the initial values to X in = (2, 2, 2), one slow cycle involving unique stick-stip motion is performed. When reseting X in = (0, 0, 0), the quasi-periodic oscillation can be observed, which can be attributed to the fact that the two subcritical Hopf bifurcations HB 2± are unreachable in this case since they occur at w = ±0.3627 (see Figure 5d1,d2).
Furthermore, the two persistence bifurcations PB ± on the two boundary equilibrium points will be unreachable when W ≤ 0.1. For instance, in the last case with W = 0.1, the slow cycle in Figure 5c2 degenerates to the unique sticking motion where the full system only behaves in a fixed point at the pseudo-equilibrium point.  Overlaps of Poincare mappings with τ = 2 k π + π/2(k ∈ N * ) and the corresponding limit cycles on the Poincare section.

Influence of the Excitation Frequency
Based on the above numerical simulations with changing excitation amplitude W, it can be found that the attractors and their bifurcations are the two key elements needed to understand the full system responses; in other word, it seems that the full system responses are entirely determined by whether the bifurcations are reachable or not. However, some interesting and different results may be observed in the full system by changing the excitation frequency Ω. In order to account for that, here we take W = 0.655 as an example to exhibit more details; one quasi-periodic oscillation in Figure 5b is on the stage with Ω = 0.001. Figure 8 gives the full system responses in the sense of the transformed phase portraits, respectively corresponding to Ω = 0.00514, Ω = 0.00515 and Ω = 0.00516, where one may find that four interesting period-1 solutions are observed with slight changes in the excitation frequency on the slow timescale.
Firstly, when Ω = 0.00514, two alternations between the stable slow manifold M 2S and the spiking attractor LC S appear in one period, leading to the fact that the full system then performs one symmetrical periodic bursting oscillation consisting of two spiking states and two quiescent states involving the unique stick-slip motions on the slow timescale, as shown in Figure 8a. Secondly, when Ω = 0.00515, only one alternation between M 2S and LC S appears in one period, leading to the two coexisting periodic bursting oscillations, respectively presented in Figure 8b with initial values X in = (0, 0, 0) and Figure 8c with initial values X in = (2, 2, 2), and only one spiking state and one quiescent state involving the unique stick-slip motion on slow timescale can be observed. In the last case, when Ω = 0.00516, one symmetrical periodic oscillation only evolving on the spiking attractor is on the stage, which certainly is not a bursting oscillation since there is no alternation between the slow manifold and the spiking attractor in the movement, as shown in Figure 8d.
As we mentioned earlier in Figure 3a, there is actually no bifurcation leading to the transition between spiking attractor LC S and stable M 2S ; moreover, note that the excitation frequency indeed does not influence the attractors as well as their bifurcations. That is to say, the transition mechanism of the periodic bursting oscillations in Figure 8a-c cannot be explained by the generalized slow fast analysis method. In addition, how to name and classify these bursting oscillations is a difficult topic since the widely applied classification method is to use the two important bifurcations leading to the transitions between spiking states and quiescent states.

Conclusions and Discussions
By introducing a low frequency excitation to a modified Chua's circuit system with a discontinuous nonlinear resistance, we built a Filippov slow-fast dynamic system. Meanwhile, the attractors as well as their conventional and nonconventional bifurcations of the fast subsystem have been analyzed by adopting property parameters based on the generalized slow-fast analysis method. It is found that the excitation amplitude and frequency are two important parameters of the system that may play an important role in the dynamical behaviors of the full system.
Firstly, taking the excitation amplitude as the analysis parameter, four different visiting modes of slow external excitation with a small enough excitation frequency are given according to the stability structures induced by conventional bifurcations. Various full system responses in every visiting mode, such as conventional sub-Hopf/LPC periodic bursting oscillation, quasi-periodic oscillations on spiking attractor, slow cycle involving unique stick-slip motion and sticking behavior, are presented in numerical simulations.
Secondly, taking the excitation frequency as the analysis parameter, although it cannot change the dynamics of the fast subsystem, the quasi-periodic oscillation in Figure 5b may degenerate to interesting dynamical evolutions from symmetrical periodic bursting oscillations to two coexisting asymmetrical bursting oscillations and finally to periodic movement on the spiking attractor via just a slight change of the excitation frequency within the slow timescale; see in Figure 8. In particular, it should be pointed out that the generation mechanism of the bursting phenomena in this case cannot be revealed by the bifurcations of the fast subsystem based on slow-fast decomposition, since there actually is no bifurcation leading to the transition between the spiking attractor and slow manifold.
Furthermore, our works show that although the unconventional bifurcations such as sliding bifurcations and persistence, which are unique in Filippov systems, cannot bring about the transition behavior between fast and slow timescales, they can affect the structure of bursting oscillatory attractors, such as the variable nonsmooth structures in the spiking state and the unique quiescent state involving unique stick-slip motion. Therefore, only considering the bifurcations leading to the transitions between spiking states and quiescent states to classify the bursting oscillations in this paper cannot be applicable here. Therefore, we may say that the insight into slow-fast dynamics in Filippov slow-fast systems is still an open problem worthy of further discussion.
In previous works [27,32], the change in the excitation frequency can affect the sparseness of the spikes in the spiking states. However, although the change in the excitation frequency cannot change the stabilities and bifurcations of the attractors in the fast subsystem, our work shows that the excitation frequency is a very important parameter. In our work, we find that a change in the excitation frequency may lead to the birth of inexplicable bursting phenomena under the special bistability structure in this paper. How to explain the generation mechanism of such bursting phenomena will be our main concern in future.