On occurrence of mixed-torus bursting oscillations induced by non-smoothness

The paper devotes to the slow–fast behaviors of a higher-dimensional non-smooth system with the coupling of two scales. Some novel bursting attractors and interesting phenomena are presented, especially the so-called mixed-torus bursting oscillations, in which the trajectory moves along different tori in turn, though the bursting attractor still behaves in quasi-periodic form. Based on a 4-D hyper-chaotic model with two scrolls, a modified non-smooth slow–fast version is established. Upon the smooth and non-smooth bifurcations analysis, the coexisted attractors with the variation of the slow-varying parameter are derived. With the increase in the exciting amplitude, bursting oscillations may change from periodic to quasi-periodic attractor, the mechanism of which is obtained by employing the overlap of the transformed phase portrait and coexisted attractors as well as the bifurcations. Sliding along the boundary on the trajectory can be observed on the bursting attractor, the mechanism of which can be accounted for by employing the non-smooth theory or by the in-turn influence of two pseudo-attractors in different regions. Super-cone defined by an unstable focus may lead to the period-adding or period-decreasing phenomena in the spiking oscillations, while fold limit cycle bifurcation may result in the dramatical change of the spiking amplitude. Different limit cycles may involve the full vector field, resulting in the so-called mixed-torus bursting oscillations, in which the quiescent states seem to be pivots of the oscillations.


Introduction
Bursting oscillations, characterized by the combination of large-amplitude oscillations (spiking state, SP) and small-amplitude oscillations (quiescent state, QS), can often be observed in the slow-fast nonlinear vector fields [1]. Generally, the bursting attractors can be classified by movement, such as periodic, quasi-periodic and chaotic bursting oscillations, or by the properties of the QSs and the SPs, such as point-cycle, cycle-cycle bursting oscillations. More commonly, they can be classified by the bifurcations at the alternations between QSs and SPs, such as fold/fold, fold/Hopf bursting oscillations [2].
The investigation of the dynamics of bursting oscillations may date back to the early twentieth century for the relaxation oscillations found in singular van der Pol equation [3]. However, it was until the Nobel Prize winners, Hodgkin and Huxley, presented a slowfast system, called Hodgkin-Huxley model, which can exhibit the activities of neuronal bursting phenomenon [4], researchers realized that the coupling of different scales is the reason for the special oscillations, which caused the behaviors of the slow-fast vector field to become one of the hot topics in nonlinear dynamics [5].
At the first stage, because of no valid theoretical method, most of the results focused on the phenomenon report of the bursting oscillations in different areas [6], or the numerical simulations for different patterns of the oscillations [7]. To reveal the mechanism of the bursting oscillations, upon the invariant manifolds [8], geometric singular perturbation theory has been developed to cope with the interaction mechanism between different scales with geometric view point [9]. For a typical autonomous slow-fast system, expressed aṡ x = f (x, y, μ), (Fast subsystem) y = εg(x, y, μ), (Slow subsystem) (1) where x ∈ R M , y ∈ R N , μ ∈ R K , while 0 < ε << 1 describes the ratio between the fast and slow scales, when ε → 0, the trajectories of (1) converge during the fast epochs to the solution of the fast subsystem [10] x = f (x, y, μ),ẏ = 0.
During slow epochs, on the other hand, the trajectories of (1) converge to the solution of 0 = f (x, y, μ),ẏ = εg(x, y, μ), which is differential-algebraic equation called the slow flow or the reduced system [11]. Upon the analysis of the behaviors of the fast and slow subsystems, geometric singular perturbation theory can qualitatively describe the mechanism of the dynamics of the full system [12]. Based on the theory, a slow-fast analysis method was proposed by Rinzel [13], in which the slow state variables are regarded as slow-varying parameters on the fact that the fast subsystem determines the quiescence and spiking oscillations, while the slow subsystem contributes to the modulation [14]. The attractors as well as the bifurcations of the fast subsystem with the variation of the slow-varying parameters can be derived, which can be used to account for the forms of quiescent states (QSs) and spiking states (SPs) as well as the transition mechanism between the two states [15]. By employing the method, a lot of bursting attractors have been reported [16], which can be further classified by the geometrical structures or the bifurcations [17].
The slow-fast analysis method can only be valid for the autonomous system with the coupling of two scales in time domain, in which obvious slow and fast subsystems can be obtained [18], while it cannot be used for the mechanism of the bursting oscillations in the non-autonomous system, for example, the periodically excited oscillators [19]. When there exists an order gap between the exciting frequency and the natural frequency, implying the coupling of two scales in frequency domain exists, slow-fast behaviors may also appear [20].
To investigate the mechanism of the phenomena, we proposed a modified slow-fast method [21]. For a typical vector field with excitation, written aṡ where x ∈ R M , μ ∈ R K , while 0 < ε << 1, we refer to the geometric singular perturbation theory and rewrite (4) in the generalized phase space aṡ x = f (x, μ, w), (Fast subsystem) w = −(εΩ) 2 w, (Slow subsystem) (5) in which w can be regarded as a generalized state variable, or in the forṁ x = f (x, μ, w), (Fast subsystem) w = A sin(εΩt), (Slow-varying parameter) (6) with slow-varying parameter w. The slow-varying parameter in the modified method is expressed in the form of an algebraic equation instead of a differential equation in the traditional method, while the fast subsystem is in generalized autonomous form [22]. By regarding the slow-varying parameter as a generalized state variable, the conception of transformed phase portrait, defined by Γ G ≡ {[x(t), w(t)] | t ∈ R} ≡ {[x(t), A sin(εΩt)] | t ∈ R}, is introduced, which can be used to describe the relationship between the state variables of the fast subsystem and the slow-varying w [23]. Combining the bifurcations of the attractors of the fast generalized autonomous subsystem with the variation of the slow-varying parameter w, one may obtain the mechanism of the bursting oscillations [24]. Based on the modified approach, several new types of bursting oscillations in periodically excited vector fields were reported, such as fold/Hopf bursting [25]. Table 1 presents the difference between the traditional slow-fast method and the modified approach.
Up to now, most of the bursting oscillations are obtained in the low-dimensional smooth systems [26], in which only codimension-1 bifurcations, such as Table 1 The difference between the traditional slow-fast method and our modified approach

Traditional method
Modified approacḣ fold and Hopf bifurcations, occur at the alternations between QSs and SPs [27]. We wonder what kinds of bursting oscillations may appear in non-smooth vector field, since not only the conventional bifurcations, but also the non-conventional bifurcations, such as the multiple-crossing bifurcations, the sliding and glazing bifurcations, may take place [28]. For a non-smooth dynamical system, the nonsmooth boundary divides the whole phase space into a set of regions, in which the trajectory is governed by different subsystems. According to the behaviors of the trajectory on the boundary, the switching boundary is partitioned into three different regions, i.e., the crossing region, the sliding region and the escaping region [29]. Non-conventional bifurcations at the boundary may result in special phenomena, most of which cannot be observed in smooth vector fields [30]. To reveal the characteristics of these non-conventional bifurcations, a few approaches have been proposed, including the differential inclusion theory [31] and the Filippov theory [32]. These theories can be employed to investigate the behaviors of the bifurcations on the boundary, such as the bifurcations related to boundary equilibrium point and the limit cycle passing across the boundary [33], which have been widely used in many low-dimensional non-smooth oscillators, especially in the planar vector fields because of the relatively simple structures of the trajectory [34]. However, it still remains a challenge how to explore the non-smooth dynamics in higher, even four, dimensional system, since there may coexist several stable attractors which may undergo different conventional and nonconventional bifurcations. Furthermore, in most published results, to simplify the analysis of the behaviors, symmetric non-smooth boundary introduced in the related papers does not violate the symmetric struc-tures of the full vector field. For example, the symmetric two-scroll chaos remains when a symmetric superplane boundary at X = 0 is introduced [35].
The main purpose of the paper aims at the characteristics of bursting oscillations in a higher-dimensional non-smooth oscillator with the coupling of two scales. We focus on the influence of the coexisted attractors and their related conventional and non-conventional bifurcations on the dynamical evolution when an asymmetric boundary introduced, which destroys the symmetry of the whole vector field.
In the following, based on a super-chaotic system, we establish a slow-fast non-smooth oscillator by introducing a slow-varying periodic excitation as well as a non-smooth boundary. Upon the analysis of the dynamical evolution with different distributions of the coexisted attractors and the bifurcations, some novel bursting attractors as well as the mechanism are presented, in which a few interesting phenomena such as the periodincreasing (decreasing) and mixed-torus structures can be found.

Mathematical model
The first hyper-chaotic attractor was reported by Rössler [36], which has more than one positive Lyapunov exponent, implying dynamics of the attractor expends in more than one direction. Since then, a lot of such oscillators have been established, among which a new 4-D vector field with two scrolls was proposed in the form [35] with w = 0 and g(x) = 0. To find the natural frequency of (7) we set the external excitation w = 0. The natural frequency can be approximated by the imaginary part of the pair of complex conjugate eigenvalues of the Jacobian matrix at the origin. Numerical simulations as well as electronic circuit experiment show that the multi-stability leads to the two-scroll hyper-chaotic attractor with the variation of the parameters. We refer to the system for the reasons as -It may display two-scroll hyper-chaotic attractor with two positive Lyapunov exponents for certain parameter values. -It can be easily implemented by electronic circuit.
-The dimension of the vector field is relatively high.
-The vector field is symmetric under the transforma- -Coexisted stable attractors can be observed with the variation of parameters.
In order to investigate the influence of the multistability and possible regular and non-regular bifurcations the dynamics when the coupling of two scales and non-smoothness exist simultaneously, here we introduce a slow-varying external excitation and a nonsmooth term, expressed in the forms w = A sin(εΩt) and g(x) = −3sgn(x − 1), respectively. Here the non-smooth function g(x) are introduced to model the influence of the automatic switch controlled by threshold in the circuit. When the exciting frequency εΩ is far less than the natural frequency of the system, denoted by ω N , the coupling of two scales in frequency domain exists. Furthermore, the asymmetric nonsmooth boundary, defined as , in which the trajectory is governed by two different subsystems, respectively.

Remark
The natural frequency of a nonlinear system may vary with its behaviors. For example, when the trajectory tends to a stable focus, it can be approximated by the imaginary part of the pair of conjugate eigenvalues, while when the trajectory oscillates around a limit cycle, it can be computed at the frequency of the cycle [37].
As we have mentioned in the above section and our other papers [38][39][40], for 0 < εΩ 1, the whole exciting term w = A sin(εΩt) can be regarded as a generalized state variable, which forms the slow subsystem, leading to the so-called fast subsystem in generalized autonomous form. In order to reveal the mechanism of the dynamics of full system, now we turn to the attractors and their bifurcations of the fast subsystem with the variation of the slow-varying parameter w.

Bifurcation of the generalized autonomous fast subsystem
The generalized autonomous fast subsystem can be further expressed in the forṁ which implies the trajectory of the fast subsystem is governed by different subsystems when it is located in different regions D ± . Therefore, we need to investigate the bifurcations of the two subsystems defined in the two regions D ± as well as the whole fast subsystem, respectively. Here for simplicity we denote the full nonsmooth generalized autonomous fast subsystem in (7) by FS 0 and its two subsystems by FS ± , located in D ± , respectively. Symmetric structures can be found in both the two subsystems FS ± since they keep unchanged under the transformation All the parameters in (7) not only affect the number and the properties of the equilibrium points, but also the bifurcation details. It is very tough to analyze the general case. To ensure coexistence of multiple equilibrium branches and their bifurcations, such as fold and Hopf bifurcations, for the fast subsystem in two different regions, here we consider a typical case with the parameters fixed at to derive the equilibrium branches and their bifurcations of the three subsystems FS 0 and FS ± with the variation of the slow-varying parameter w.

Equilibrium state and the bifurcations of the subsystems FS ± in D ±
For the subsystem in D + with g(x) = −3, the equilibrium point can be written as where X 1+ satisfies while X 2+ = + , with Similarly, the equilibrium point in D − can be expressed as while The stability of the equilibrium points can be uniformly determined by the associated characteristic equation, written in the form where The equilibrium point is stable if Re(λ i ) < 0, i = 1, 2, 3, 4, which may lose its stability via two types of codimension-1 bifurcations. Fold bifurcation may occur when one zero eigenvalues appears, resulting in the jumping phenomenon between different equilibrium points, while Hopf bifurcation may take place when a pair of pure imaginary eigenvalues can occur, leading to the periodic oscillations. Because of the complicated expressions of the equilibrium points as well as the characteristic equation, we now turn to the numerical simulations. Figure 1 gives the steady states of the generalized autonomous fast subsystem FS 0 and the two subsystems FS ± in D ± with the variation of the slow-varying parameter w.
We first focus on the subsystem FS + in D + . With the variation of w, the equilibrium points form three equilibrium branches, denoted by EB + i , i = 1, 2, 3, see Fig. 1. Because of the symmetry of FS + , Table 2 only presents the dynamical behaviors for w ≤ 0, while the details for w > 0 is omitted here for simplicity.
Note All the equilibrium points on EB + 1 are pseudoequilibrium points, while the number of equilibrium point on E + 3 varies from one to three, which may appear as admissible or pseudo-equilibrium points, depending on the relationship between their locations and the nonsmooth boundary.
For the subsystem FS + in D + , a stable node, an unstable focus and an unstable saddle can be observed Two stable node and a saddle EB + when w ∈ (−∞, −2.1218), shown in Fig. 2 for w = −6.82. Two hetero-clinic orbits, denoted by S 1,2 , can be observed. S 1 connects the unstable focus E + 1 and the stable node E + 3 , while S 2 connects the saddle point E + 2 and the stable node E + 3 . Computation reveals that at w = −2.1218 and w = −1.2394, subcritical Hopf bifurcations occur, causing the qualitative changes of the phase portrait on the phase space. When w ∈ (−2.1218, −1.2394), a stable focus E + 1 , an unstable saddle E + 2 and a stable node E + 3 can be observed, where two hetero-clinic orbits S 1,2 connect E + 2 and E + 3 , shown in Fig. 3 for w = −1.50.
When w ∈ (−1.2394, −1.2386), the stable focus E + 1 changes to an unstable one via subcritical Hopf bifurcation, which meets the saddle point E + 2 at w = −1.2386 to form a saddle-focus bifurcation, leading to the disappearance of the two equilibrium points.
When w ∈ (−1.2386, −0.1560), only one stable node E + 3 exists. Fold bifurcation in type of saddle-node at w = −0.1560 results in the occurrence of a saddle Symmetrically behaviors can be observed for w ∈ (0, +∞). For example, subcritical Hopf bifurcations at w = +1.2394 and w = +2.1218 lead to the alternation of E + 1 between a stable and an unstable focus, see Fig. 5 for w = +1.50, with two hetero-clinic orbits S 1,2 connecting E + 2 and node E + 3 , respectively.
Remark Though the dynamical behaviors are symmetric about w = 0.0, some of the equilibrium points may change between pseudo-and admissible equilibrium points because of the locations of the points may change with the variation of w.
Now we turn to the subsystem FS − in D − . With the variation of w, the equilibrium points form also three equilibrium branches, denoted by EB − i , i = 1, 2, 3, see  Table 3 presents the equilibrium points as well as the bifurcations of FS − with the variation of w for w ≤ 0, while the details for w ≥ 0 are omitted because of the symmetry. For w ∈ (−∞, −4.4363), a stable node and an unstable node on EB − 2 can be observed, which is coexisted with a stable limit cycle L 1 bifurcated from an unstable focus on EB − 1 , see Fig. 6a for w = −6.82. The limit cycle may disappear via super-critical Hopf bifurcation at w = −4.4363, at which the unsta-ble focus changes to be a stable one, see Fig. 6b for w = −3.0. It can be observed that there exist two hetero-clinic orbits, denoted by S 1,2 , where S 1 connects the stable focus E + 1 and the unstable node E + 2 , while S 2 connects the point E + 2 and the stable node E + 3 . With the increase in w, E + 2 moves toward E + 3 , and finally meet E + 3 at w = −2.9179 to form a nodenode bifurcation, resulting in the disappearance of the two nodes, which leads to only a stable focus for w ∈ (−2.9179, 0.0).
Remark Some of the equilibrium points may change between pseudo-and admissible equilibrium points because of the locations of the points may change between D ± with the variation of w. Furthermore, because of the high dimension of the system, it is very difficult to obtain exactly the attracting basins of the stable attractors.
3.2 Bifurcation of the full non-smooth generalized autonomous fast subsystem FS 0 Now we turn to the full non-smooth generalized autonomous fast subsystem FS 0 . From the non-smooth boundary h(x) = x = 1, we can compute the directional Lie derivative as L F ± =< ∇h, F ± >= f ± 1 . Setting L F ± = 0 yields the set of tangent singular points, Therefore, for a given value of w, when the trajectory reaches the non-smooth boundary at the point X 0 ∈ Σ, if f ± 1 = 0, sliding or grazing phenomenon can be observed. Otherwise, the trajectory may pass across the boundary. Furthermore, with a small perturbation of the parameters, if the value of f ± 1 = 0 keeps at zero, the trajectory may sliding along the boundary; otherwise, a grazing point on the trajectory can be observed. Now we turn to the numerical simulations for nonsmooth bifurcations. Table 4 shows the attractors in the FS 0 with the variation of w. Figure 7 gives the solutions of FS 0 and FS ± for w = −7.0. For the subsystem FS + in D + , there exists a hetero-clinic orbit S connecting the unstable focus E + 1 and the stable node E + 3 . While for the subsystem FS − in D 2 , the trajectory may approach a stable limit cycle L 1 , which have two interaction points M 2 and M 3 with the non-smooth boundary Σ. When the trajectory of FS 0 is located in D − , it is governed by FS − and tends w ∈ (−∞, −6.8203) A stable limit cycle L 3 and a stable node on EB + 2 w = −6.8203 C-Bifurcation w ∈ (−6.8203, −4.5433) Two stable limit cycles L 1 an L 3 , together with a stable node on EB + 2 and an unstable limit cycle L 2 w = −4.5433 Fold limit cycle bifurcation of L 2 and L 3 w ∈ (−4.5433, −4.4363) A stable limit cycle L 2 and a stable node on EB + Grazing bifurcation w ∈ (+1.6990, +3.0656) A stable grazing limit cycle L 5 w = +3.0656 Period-doubling bifurcation w ∈ (+3.0656, +4.7195) Period-doubling oscillations w = +4.7195 Period-doubling bifurcation w ∈ (+4.7195, +∞) A stable crossing limit cycle L 5 to the stable limit cycle L 1 . When it arrives at the point M 2 in Fig. 7a, since dx 1 /dt keeps positive, see Fig. 7b, it passes across the boundary to enter D + . Then, the trajectory turns to be governed by FS + . Affected by the hetero-clinic orbit, the trajectory moves along the orbit of the unstable focus. When it returns to boundary, it also passes across because of the same reason. The alternated influence may finally reach a balance, leading to a non-smooth limit cycle L 3 in FS 0 . Figure 8 gives the phase portraits on the (x 1 , x 2 ) plane of FS 0 for w = −5.0 and w = +2.0, respectively. Note that when w = −5, 0, the stable limit cycles L 1 and L 3 coexist, corresponding to different attracting basins. However, because of the grazing From the enlarged phase portrait in Fig. 9a, one may find that there exist two types of the behaviors when the trajectory arrives at the non-smooth boundary. It may graze the non-smooth boundary, seeing point M 1 in Fig. 9a. For the case, one may find from the time history of x 1 in Fig. 9b, c that the value of x 1 at M 1 reaches its extreme value, implying dx 1 /dt = 0, which agrees well with the theoretical result. It may also pass directly across the boundary, seeing points M 2 and M 3 in Fig. 9a, at which the related values of dx 1 /dt keep the sign, demonstrated by the time history of x 1 in Fig. 9c.
With the increase in the parameter w, perioddoubling bifurcation can be observed for FS 0 , see  Further increase in w may cause the disappearance of the grazing point, leading to period-1 oscillations, shown in Fig. 11 for w = +6.0. Period-adding phenomenon Numerical simulation reveals that both the periods of the two limit cycles L 1 and L 3 change slowly with the variation of w, which is almost similar to the phenomenon observed for the limit cycle in smooth vector field. However, the period of the limit cycle L 5 changes dramatically with the variation of w. The limit cycle L 1 via Hopf bifurcation at H B − 1 in Fig. 1 is governed by the subsystem FS − in D − , the amplitude of which increases with the decrease in w. When it reaches the boundary, C-bifurcation occurs, leading to the stable non-smooth limit cycle L 3 , which causes the abrupt increase in the oscillating amplitude. However, small change of the period can be observed between the two cycles L 1 and L 3 , since most part of the non-smooth limit cycle L 3 is still governed by the subsystem FS − , while the other part is determined by the flow around the focus E + 1 in D + , see Figs. 2b and 7. Now we focus on the limit cycle L 5 via grazing bifurcation at the boundary. From the stability analysis of the equilibrium branches, one may find, when the trajectory is located in region D − , it may passes across the boundary to enter region D + . Once the trajectory is located in D + , it will be governed by the subsystem FS + , the phase portrait of which is symmetric to that in Fig. 2. Before the period-doubling bifurcation occurs, the trajectory is mainly governed by the subsystem FS + except in the small neighborhood of the grazing point. Note that the flow surrounding the unstable focus on EB + Fig. 12 Comparison between the three cases. a Phase portrait on the (x 1 , x 2 ) plane; b Time history of x 1 adding/decreasing phenomenon with the slow-varying w (Fig. 12). Period-doubling phenomenon For the value of w, there exist two hetero-clinic orbits connecting a stable in D − and an unstable node in D + , which implies there exist two possible routes for the trajectory from any initial point located in D + to approach the stable node in D − , when the non-smooth boundary does not exist, shown in Fig. 5. Therefore, two possible intersection points between the trajectory and Σ may be observed. When only one of the possible intersection points appears, which behave in grazing type, period-1 solution can be observed, seeing N 3 in Fig. 13a for w = +2.0, while the other point just approaches the boundary, seeing N 1 in Fig. 13a. With the increase in w, the inertia of the movement increases, causing the original N 1 to arrive at the boundary. Because of the characteristics of the vector field, no singular property of the intersection point exists, which causes the trajectory to pass directly across the boundary, yielding two more intersection points on the boundary, resulting in period-doubling bifurcation, see Fig. 13b for w = +4.0.
Further increase in w may continue to increase the vibration intensity, causing the trajectory at N 3 to directly pass across the boundary, which lead to the disappearance of the grazing point on the trajectory. For the situation, period-1 solution can be observed, which implies the trajectory may pass across the boundary only via the two possible points.
Remark Here we would like to point out that the supercone determined by the stable or unstable focus may be one of the reasons for the period-decreasing or periodincreasing phenomenon of the periodic oscillations in the non-smooth vector field. Furthermore, the grazing point may be eliminated by the increase the vibration intensity, which can be used for controlling purpose.
From the above analysis, one may find that for a non-smooth vector field, there exist different attractors as well as regular and non-smooth bifurcations, which leads to different distributions of the attractors and the bifurcations with the variation of w. In the following, we will investigate the slow-fast behaviors when the coupling of two scales in frequency involves the full system.

Evolution of the dynamics
Note that the slow-varying parameter w may change between −A and +A periodically for a given exciting amplitude A. With the increase in the exciting amplitude, different types of attractors and their bifurcations may involve the full system, resulting in qualitatively different behaviors. Now we focus on the dynamical evolution with the increase in the exciting amplitude.
Since the distance between some bifurcation points is very short, to display the influence of all possible stable attractors of the generalized autonomous fast subsystem on the dynamics of the full system, we need the variable w to change slowly enough. For the purpose, here we fix the exciting frequency εΩ at εΩ = 0.001.
Remark Numerical simulations demonstrate that when the small parameter εΩ is taken a bit large, some of the behaviors may disappear, since w passes across the region confined by the associated bifurcation points with short distance along w-axis before the influence of the bifurcations occurs. It can be found that the stable EB − 1 meet the nonsmooth boundary at w = 1.5536, see Fig. 1. Therefore, for A ∈ (0, 1.5536), though there exist a few bifurcations, such as node-node bifurcation and the subcritical Hopf bifurcation, the full system appears in limit cycle with the same frequency as that of the excitation, since the trajectory of the full system does not pass across the boundary Σ, implying the full system oscillates around the equilibrium branch EB − 1 governed by the subsystem in D − , which we omit here for simplicity, see Fig. 1.

Periodic relaxation oscillations with sliding
For A = 1.5536, EB − 1 reaches the non-smooth boundary, causing the trajectory to slide along the non-smooth boundary, see Fig. 14 for A = 1.60.
The dynamics behaves in periodic oscillations with the same period as that of the excitation, the trajectory for each period of which can be divided into three segments, denoted by P 1 P 2 , P 2 P 3 and P 3 P 4 , where the segment P 2 P 3 is located on the boundary while repetitive spiking oscillations can be observed around the point P 4 .
To reveal the mechanism of the non-smooth bursting oscillations, we turn to the overlap of the transformed phase portrait and equilibrium branches as well as their bifurcations on the (w, x 1 ) plane, shown in Fig. 15a.
With the increase in the slow-varying parameter w, the trajectory, starting from the point P 1 with minimum value of w at w = −1.60, moves almost strictly along the stable segment of equilibrium branch EB − 1 , behaving in quiescent state, since the trajectory is governed by the subsystem in D − . When the trajectory reaches the non-smooth boundary Σ at P 5 , it slides along Σ because of the sliding bifurcation. When w increases to the maximum value with w = +1.60 at P 2 , the tra- jectory turns to the left. It then moves almost strictly along Σ to the point P 3 , at which non-smooth sliding bifurcation occurs. The trajectory jumps to the stable segment of equilibrium branch EB − 2 , yielding repetitive spiking oscillations. The trajectory finally settles down to the stable equilibrium branch at the point P 4 , and turns to move almost strictly along the stable equilibrium branch. When the trajectory arrives at the point P 1 , one period of the oscillations is finished.
The structure of the bursting oscillations is presented in Fig. 15b, which can be can be called periodic nonsmooth relaxation oscillations with sliding, because of the sliding bifurcation at the transitions between the boundary and the stable equilibrium branches. Spiking oscillations When the trajectory along the boundary enters the region D − , it may tend to the stable equilibrium branch of the subsystem in D − , which is of focus type. Therefore, the spiking oscillations can be approximated by the transient process from the point on the boundary to the focus on the equilibrium branch EB − 1 for the corresponding value of w. Sliding mechanism Once the trajectory moving along the stable equilibrium segment EB − 1 of the subsystem in D − reaches the non-smooth boundary Σ at P 5 , it then turns to move along the boundary, resulting in sliding phenomenon. From the time history of x 1 and w in Fig. 14b, one may find when the trajectory along EB − 1 reaches the boundary at P 5 with x 1 = 1.0, the value of w at P 5 can be approximated at w = 1.5818, which implies w may still increase with the evolution of time since A = 1.60. Theoretically, without the boundary, the trajectory passes across the boundary and enters the region D + . Because of the boundary, when the trajectory passes across the boundary, it turns to be governed by the subsystem FS + in D + . Since the associated stable equilibrium branch EB + 2 of the system in D + is Fig. 16 Sliding mechanism located in D − , though it is a pseudo-stable equilibrium branch, it still attracts the trajectory, causing the trajectory to return to the region D − , see Fig. 16. When the trajectory passes across the boundary to enter the region D − , it then is governed by the subsystem FS − in D − . Note that for w ∈ (1.5818, 1.60], the stable equilibrium points EB − 1 of the subsystem FS − is located in D + , which causes the trajectory to return to region D + . Skipping the transient process, the back and forth influence of the two governing subsystems FS ± finally causes the trajectory to move along the non-smooth boundary, leading to the sliding phenomenon. When the slow-varying parameter w decreases from the maximum value w = 1.60, because of the inertia of the movement along the boundary, the trajectory may still move along the boundary for a segment until the effect of sliding bifurcation occurs, see Fig. 16. The sliding mechanism can also be demonstrated by the non-smooth theory. When w ∈ (−1.5818, 1.60), for the trajectory in D + , x 1 may decrease because of the stable equilibrium branches EB + 2 , while for the trajectory in D − , x 1 may increase because of the stable EB − 1 . Assuming the trajectory from D ± + arrives at Σ at the point X ± , one may find f + 1 (X + ) < 0 and f − 1 (X − ) > 0, which implies the segment w ∈ (−1.5818, 1.60) on w-axis corresponds to the sliding region, with sliding bifurcation at P 5 . sliding bifurcation at P 5 because of the slow passage influence caused by the inertia of the movement along the boundary, which causes the jumping phenomenon to occur at P 3 instead of P 5 .

Periodic point-cycle bursting oscillations
With the increase in the exciting amplitude, when A = 1.6990, grazing bifurcation at w = +1.6990 may involve the dynamics, leading to the qualitative change of the behavior, see Fig. 17 for A = 1.70. From the phase portrait in Fig. 17a, one may find that the topological structure of the part of the trajectory located in the region D − is the same as that in Fig. 15a, since the trajectory is governed by the same subsystem FS − . However, instead of sliding along the boundary, for larger exciting amplitude, the trajectory, not only slides along the boundary, but also passes across the boundary to oscillate according to a cycle, yielding repetitive spiking oscillations, seeing SP in Fig. 17c, since stable non-smooth limit cycle meets the boundary at w = +1.6916 to yield grazing bifurcation.
From the time history in Fig. 17b, one may find that, the trajectory for one period can be divided into four segments, corresponding to two quiescent state and two spiking state. The trajectory in spiking state for the part P 3 P 4 , denoted by SP 1 , oscillates according to the non-smooth limit cycle, while the trajectory in the part (a) (b) Fig. 18 Overlap of the transformed phase portrait and equilibrium branches on the (w, x 1 ) plane in a with locally enlarged part in b for A = 1.70 P 6 P 5 , denoted by SP 2 , corresponds to the transient process from the point on the boundary to a stable focus. The quiescent QS 1 for the trajectory from P 1 to P 3 is composed of two parts, where the trajectory in P 1 P 2 moves along the stable equilibrium branches of the subsystem in D − , while in P 2 P 3 , it slides along the boundary. The trajectory in quiescent state QS 2 from P 4 to P 6 just moves along the boundary, shown in Fig. 17c, d. Furthermore, the frequency of SP 1 can be approximated by the frequency of the non-smooth limit cycle at Ω 1 = 0.5148, while the frequency of SP 2 can be approached by the imaginary part of the pair of complex conjugate eigenvalues at Ω 2 = 2.0944, both of which agree with the numerical results. For example, from the time history in Fig. 17d, the frequency of SP 1 can be computed at Ω 1 = 2π/T 3 = 0.5236.
We also turn to the overlap of the transformed phase portrait and equilibrium branches on the (w, x 1 ) plane to investigate the mechanism of the bursting oscillations, see Fig. 6 (Fig. 18).
When the trajectory starts from the point P 1 , at which w = −1.70, it moves almost strictly along the stable equilibrium branch EB − 1 , appearing in quiescent state QS 1 , until it reaches the boundary at P 2 , at which sliding bifurcation occurs. With the increase in w, the trajectory moves along the boundary, leading to the sliding phenomenon. When the trajectory arrives at the point P 3 , it turns to oscillate according to the limit cycle L 5 , leading to spiking state SP 1 via grazing bifurcation.
When the trajectory oscillates to P 7 with maximum w at w = +1.70, it turns to the left in spiking oscillations since w decreases with the increase in time. Grazing bifurcation at P 4 causes the trajectory to settle down to the non-smooth boundary to begin the quiescent state QS 2 . When the trajectory moves along the boundary to the point P 6 , sliding bifurcation occurs, causing the tra- jectory to jump to the stable equilibrium branch EB − 1 in D − , yielding spiking oscillations SP 2 . The trajectory may finally settle down to the stable equilibrium branch at P 5 and turns to move almost strictly along the branch EB − 1 until the trajectory arrives at the starting point P 1 , which finishes one period of the bursting oscillations.
The structure of the bursting oscillations is shown in Fig. 19, from which one may find that sliding bifurcation and grazing bifurcation cause the alternations between the quiescent states and spiking states. Therefore, the type of busting oscillations can be called nonsmooth point-cycle bursting from the structure, or be called sliding/G-bifurcation bursting attractor by the bifurcations at the transitions.

Bursting attractor with one-torus spiking oscillations
More influence of the non-smooth grazing limit cycle L 5 on the bursting attractor can be found for larger exciting amplitude, see Fig. 20 for A = 2.50 and A = 3.0. The overlap of the transformed phase portrait and the equilibrium branches for A = 3.0 is presented in Fig. 21a, while the related time history of x 1 and w is shown in Fig. 22.
From the overlap of the transformed phase portrait and equilibrium branches on the (w, x 1 ) plane in Fig. 21a, one may find that the mechanism of the bursting oscillations is the same as that of the attractors in Fig. 17. The repetitive spiking state SP 2 oscillate Numerical simulations reveal that the Poincaré map behaves in a fixed point when the initial condition is taken at an arbitrary point on the trajectory not in the spiking oscillations SP 2 . However, when the initial condition is taken at a point on the trajectory in the spiking oscillations SP 2 , the projections of the Poincaré map on the (x 1 , x 2 ) plane are located exactly on a cycle, which indicates that the spiking oscillation SP 2 appears in quasi-periodic form, shown in Fig. 21b. Period-decreasing and period-adding in spiking oscillations Now we focus on the time history of the x 1 , presented in Fig. 22. When w takes at its minimum value at P 1 with w = −3.0 in Fig. 22a, the trajectory moves almost strictly along EB − 1 to the non-smooth boundary. Then, the trajectory turns to slide along the boundary via sliding bifurcation, until grazing bifurcation occurs, which causes the trajectory to begin repetitive spiking oscillations. For the spiking oscillations SP 1 , it can be found that the slow-varying parameter w may change between w = +1.6990 to w = +3.0, which implies that spiking state may oscillate according to different non-smooth limit cycles. Though the amplitudes may vary almost strictly according to L 5 , the frequencies change with the variation of w, shown in Fig. 22.
The period of the spiking oscillations may decrease when w increase from w = 1.6990. For example, when w = 2.1040, the period can be approximated at T 1 = 4.0513, shown in Fig. 5c, while for w = 2.9839, the period is T 2 = 2.3981, see Fig. 5d. Further investigation reveals that period-decreasing phenomenon can be observed when w increase from w = +1.6990 to w = +3.0.
On the contrary, one may find the period can be computed at T 3 = 2.3984 when w = 3.0, see Fig. 22e, while T 4 = 9.1890 for w = 1.6701, shown in Fig. 22f, which implies period-adding phenomenon can be obtained for w decreases from w = +3.0 to w = +1.6990.
When the trajectory oscillates to the boundary at P 4 in Fig. 22f, it turns to move along the boundary to start the quiescent state, yielding sliding phenomenon until the effect of the sliding bifurcation occurs, which causes the trajectory to jump to the stable focus-type equilibrium branch EB − 1 .
Remark Not only the period-decreasing and periodadding phenomena, but also the grazing forms, are similar to those of the non-smooth fast subsystem, see Figs. 8b and 12, which further demonstrate that the super-cone caused by the unstable focus is one of the reasons for the period-adding and period-decreasing phenomena in non-smooth vector field. Furthermore, when the initial condition is not taken at the SP 1 stage, the projections of Poincaré map on the (x 1 , x 2 ) plane will be located at a fixed point. Therefore, both the bursting attractors in Fig. 20 are of quasi-periodic type, in which the quiescence seems to be pivot of the oscillations, implying for each period of the external excitation, the trajectory will repeat the movement in quiescent stage, while it may visit different points on the limit cycle L 5 during the spiking stage.

Bursting attractor with two-torus spiking oscillations
Note that w ∈ [−A, +A]. For A ∈ (3.0656, 4.4363), because of the interaction between the trajectory and the non-smooth boundary, period-doubling bifurcation of the non-smooth limit cycle L 5 in FS 0 can be observed, seeing the bursting oscillations in Fig. 1, which causes the qualitative change of the bursting oscillations, shown in Fig. 23 for A = 4.0.
The spiking oscillations SP 1 in Fig. 23b can be further divided into two segments approximated by w = w pd ≈ 3.2622 , in which the part R corresponds to the oscillations according to the period-2 limit cycle bifurcated from L 5 , while the other oscillate according to the limit cycle L 5 , see Fig. 23b.
From the locally enlarged parts of time history of x 1 plotted in Fig. 24, one may find the frequency of repetitive spiking oscillations in SP 1 changes with the variation of w. When the trajectory begin to spike at P 2 in Fig. 24a, the frequency of the oscillations can be computed at ω 1 = 2π/T 1 = 0.7893. It increases with the increase in w because of period-decreasing phenomenon until period-doubling bifurcation occurs at w = w pd with the frequency approximated at 2.7312. The frequency changes to be ω 2 = 2π/T 2 = 1.3962 at w = 3.217, shown in Fig. 24b, which increases to ω 3 = 2π/T 3 = 2.02683 at w = 4.0, shown in Fig. 24c. Therefore, the period-decreasing or increasing phenomenon not only exists in the oscillations according to L 5 , but also exists in the oscillations according to the period-2 cycle bifurcated from L 5 .
The mechanism of the bursting attractor is similar to that for the A = 3.0 except that period-doubling phenomenon exists in SP 1 . From the overlap of the transformed phase portrait and the equilibrium states of the fast subsystem in Fig. 25, one may find, when the trajectory starting from P 1 in Fig. 25a  it moves almost strictly along EB − 1 to the boundary at P 2 in Fig. 25b. Sliding bifurcation causes the trajectory to move along the boundary until grazing bifurcation occurs P 3 in Fig. 25b, at which spiking oscillations SP 1 according to L 5 take place.
The trajectory oscillating according to L 5 turns to oscillate according to the period-2 cycle when perioddoubling bifurcation occurs at P 4 in Fig. 25c. When the trajectory returns when w reaches its maximum value with w = +4.0, see Fig. 25d. Period-doubling bifurcation causes the trajectory to move oscillate according to L 5 until grazing bifurcation takes place at P 3 , which causes the trajectory to sliding along the boundary. When the trajectory moves to P 5 , the effect of sliding bifurcation appears, leading to the jumping phenomenon. The trajectory finally settles down to EB − 1 at P 6 in Fig. 25b. When the trajectory moves along EB − 1 to the starting point P 1 , one period of the excitation is finished.
To reveal the property of the oscillations, we also introduce the Poincaré map defined by N = 0, 1, 2, . . .}. When the initial point is taken at the point on the trajectory in quiescent state, the Poincaré map appears in a fixed point, which demonstrate that the trajectory will repeat the movement of quiescence for each period of the external excitation. Here we take two cross sections for the trajectory in spiking state SP 1 , corresponding to period-1 and period-2 limit cycles, respectively, shown in Fig. 26 (Fig. 27). a 9275, −1.0682, 0.9604, 0.8542, 13360.) 2.0113, −1.1731, 1.1562, 0.8909, 14118.). The cross section Π 1 is located at the region associated with the stable period-1 limit cycle L 5 , while the cross section Π 2 corresponding to period-2 cycle bifurcated from L 5 via period-doubling bifurcation The projections of the Poincaré maps on the cross section Π 1 and Π 2 are presented in Fig. 26, which show that the quasi-periodic spiking oscillations change between two tori in turn, divided by the period-doubling bifurcation point at w = w pd .
Remark By skipping long enough the transient process, for each period of the excitation, numerical simulations show that the trajectory visits the same point it is in quiescent state. However, the trajectory may visit different points located on the related period-1 and period-2 limit cycles, implying the spiking stage is composed of two segments, during which the trajectory oscillates according to two limit cycles. Furthermore, the trajectory in SP 2 is repeated in each period, which can be approximated by the transient process from the point P 5 on the boundary to the stable focus P 6 on EB − 1 . The trajectory oscillates around two cycles, see Fig. 28a, which can be divided into six segments corresponding to three stages of quiescent states and spiking states, respectively, denoted by QS i and SP i with i = 1, 2, 3, shown in Fig. 28b.
The mechanism can be obtained by the overlap of the transformed phase portrait and the equilibrium states of the fast subsystem, plotted in Fig. 29. The trajectory, starting from P 1 with w = −5.0, behaves in spiking oscillations SP 3 , which settles down to EB − 1 at P 2 via super-critical Hopf bifurcation to begin the quiescence QS 1 , shown in Fig. 29a.
The trajectory moves almost strictly along EB − 1 to the boundary Σ at P 3 , at which sliding bifurcation occurs, causing the trajectory to sliding along Σ to P 4 . It turns to spike according to the non-smooth limit cycle L 5 via grazing bifurcation, represented by SP 1 . The trajectory may oscillate according to a period-2 cycle dur- ing the interval between P 5 and P 6 , corresponding to the period-doubling and reverse period-doubling bifurcation points, respectively, shown in Fig. 29a.
The trajectory then oscillates according to L 5 until w reaches its maximum value at P 7 with w = +5.0. Similar behaviors can be observed for the trajectory to return from P 7 to P 4 . The trajectory changes to slide along Σ via grazing bifurcation. When the trajectory in QS 2 state moves to P 8 along Σ, it jumps to EB − 1 for the effect of sliding bifurcation.
Repetitive spiking oscillations SP 2 can be observed, which may quickly settle down to the focus-type EB − 1 to begin the quiescence QS 3 . It turns to oscillate according to L 1 via Hopf bifurcation at P 1 to begin the spiking oscillations SP 3 . When the trajectory arrives at P 1 , one period of excitation is finished, see Fig. 29b.
The structure of the bursting attractor can be described by Fig. 30, from which one may find that the trajectory may oscillate around a stable focus and two limit cycles. Though one of the cycles may change between period-1 and period-2 types, and sliding phenomenon along the boundary can also be observed, the form of bursting oscillations can be called quasiperiodic 1-point/2-cycle bursting attractor from the structure. Of course, we can also name the bursting oscillations by the bifurcations at the alternations between the quiescent states and spiking states.
Further investigation reveals that the trajectory repeats its route when it is in quiescent states. However, it does not repeat its route when the trajectory behaves in SP 1 and SP 2 , seeing the locally enlarged part of the transformed phase portrait plotted in Fig. 29b, similar to that in Fig. 25d. To reveal the property of the spiking oscillations, we chose three cross sections, denoted by Pi i (i = 1, 2, 3), shown in Fig. 28b. The projections of the Poincaré map at the three cross sections, defined by Fig. 31.
It can be found that the projections of the map at the three cross sections are located on different cycles, respectively, which demonstrates that both SP 1 and SP 3 are quasi-periodic. Furthermore, SP 1 can be divided into two segments, corresponding to the quasi-periodic oscillations associated with period-1 and period-2 cycles, respectively.
Remark The phenomenon does not violate the theorem of the existence and uniqueness of a solution of a differential equation, since the trajectory visits the same point in the phase space with different time.

Bursting oscillations with four-torus spiking oscillations
Further increase in the exciting may cause the influence of C-bifurcation at w = −6.8203 involve the bursting oscillations, which may change the bursting structure, see Fig. 32 for A = 10.0. The parts related to the original two limit cycles L 1 and L 5 for A = 5.0 develop to occupy more phase space when w = 10.0. The trajectory for A = 10.0 can also be divided into six segments, where the parts QS 1,2,3 and SP 1,2 behave in qualitative the same as those for A = 5.0, while great change can be observed between the segment SP 3 comparing the two cases for  (Fig. 33). Now we focus on the SP 3 in Fig. 32b. From the locally enlarged part of the time history presented in Fig. 32d, one may find SP 3 can be further divided into two segments by the C-bifurcation, which oscillate according to L 1 and L 3 , respectively, seeing the projections of the Poincaré may on the (x 1 , x 2 ) plane for the two segments in Fig. 34.
Further simulation reveals that the frequency in SP 3 changes very slowing with the evolution of time, which agrees well with the result related to the properties of L 1 and L 3 .
Remark The C-bifurcation results in the abrupt change of the amplitude of the spiking oscillations, which can be checked by the theory of non-smooth dynamics, referring to dynamical behaviors of the non-smooth generalized autonomous fast subsystem FS 0 . Now we employ the overlap of the transformed phase portrait and the equilibrium states of the fast subsystem, plotted in Fig. 34, to investigate the mechanism of the oscillations.
The trajectory, starting from P 1 with w = −10.0, oscillates according to L 3 in spiking state SP 3 to the point P 2 , at which C-bifurcation occurs, causing the trajectory to oscillate according to L 1 , shown in Fig. 34. Dramatic decrease in the oscillating amplitude can be observed. Further increase in w leads to the gradual decrease in the amplitude of the spiking oscillations, until Hopf bifurcation takes place at P 3 , causing the trajectory to move almost strictly along EB − 1 . Since the movement of the trajectory from P 3 is similar to that for A = 5.0 until it returns to P 3 , which we omit here for simplicity. Therefore, the structure of the attractor can be described in Fig. 35. Numerical simulations reveal that the projections of the Poincaré of the attractor on the (x 1 , x 2 ) plane may change between five forms, i.e., a fixed point, two limit cycles in L 1 and L 3 forms, and a period-1 or period-2 limit cycle associated with L 5 , depending on the initial conditions. Therefore, it is a quasi-periodic bursting attractor associated with four tori (Fig. 36).
Hopf bifurcation and sliding bifurcation, together with the grazing bifurcation yield the alternations between quiescent states and spiking states on the trajectory. Therefore, the type of movement can be called quasi-periodic Hopf/sliding/grazing bursting oscillations by the bifurcations. Note that the trajectory may oscillate around the four limit cycles L 1 , L 3 , L 5 and the cycle bifurcated from L 5 via period-doubling bifurcation, together with the equilibrium branch EB 1 1 and the non-smooth boundary Σ, it can also be called quasiperiodic non-smooth 1-point/4-cycle bursting attractor.
Further increase in the exciting amplitude cannot change the topological structure of the bursting oscillations since no further bifurcation may involve the vector field, which can also be demonstrated by the numerical simulations.

Conclusions and discussion
For a relative high-dimensional non-smooth system, there exist different coexisted attractors, which may undergo different conventional and non-conventional bifurcations. When an external excitation is introduced with the exciting frequency far less than the natural frequency, several novel bursting attractors may appear because of the influence of the coexisted attractors and the bifurcations. For the model in the paper, when the trajectory does not pass across the boundary, the system behaves in periodic oscillation with no slow-fast effect since it is governed by a smooth subsystem. With the increase in exciting amplitude, sliding bursting attractor can be observed because of the sliding bifurcation of the boundary equilibrium point, which may evolve to periodic point-cycle bursting and quasi-periodic pointcycle bursting with the further increase in the exciting amplitude. The quasi-periodic bursting may behave in mixed-torus type, in which the number of torus may change from one to four, depending on the number of stable limit cycle involved. Because of the supercone structure defined by the flow around the unstable focus, period-adding/decreasing phenomenon can be observed when the trajectory oscillates according to set of super-cones surrounding an unstable focus-type equilibrium branch. Here we would like to suggest that for a periodically excited system with the coupling of multiple scales, we need to use different initial conditions to compute the Poincaré map for the property of the bursting oscillations, since there may exist different forms of Poincaré map.