Global dynamics of a SD oscillator

In this paper, we derive the global bifurcation diagrams of a SD oscillator which exhibits both smooth and discontinuous dynamics depending on the value of a parameter a. We research all possible bifurcations of this system, including Pitchfork bifurcation, degenerate Hopf bifurcation, homoclinic bifurcation, double limit cycle bifurcation, Bautin bifurcation and Bogdanov–Takens bifurcation. Besides, we show that the system has five limit cycles, including four small limit cycles and one large limit cycle. At last, we give all numerical phase portraits to illustrate our results.


Introduction and main results
In recent years, SD (smooth and discontinuous, for short) oscillator was proposed and investigated for studying the transition from smooth to discontinuous dynamics, see, for instance [2][3][4][5]15]. In [5,15], the van der Pol damped SD oscillator is given bÿ for studying this transition, where a ≥ 0, b and ξ can take arbitrary real values. More precisely, the smooth dynamics appears when a > 0, while the discontinuous dynamic behavior occurs at a = 0. The global dynamics was completely studied in [5] when a = 0, and in [15] when |a − 1| < ε, |ξ | < ε and ε are sufficiently small. Besides, Cao et al. in [2][3][4] studied the following linear damped SD oscillator with periodic excitation using the fact that when f 0 = 0 this system admits a Hamiltonian formulation, and they analyze the dynamics of the perturbed Hamiltonian system by Melnikov methods when f 0 is small. Moreover, they gave many numerical results. Clearly, the dynamics of periodic system (2) is different from the dynamics here studied, the one of system (1). Clearly, system (1) can be rewritten as the 2dimensional differential system whereξ = ξ/3,b = 3b, and for simplicity in what follows we still denoteξ andb by ξ and b, respectively. Note that system (3) is invariant by the transformation (y, t, ξ) → (−y, −t, −ξ). Therefore, we only need to consider the set of parameters where R + = [0, +∞).
In this paper, we call large limit cycles to the ones surrounding all three equilibria and small limit cycles the ones surrounding a single equilibrium. For the notions and definitions which appear in the statement of Theorem 1, see its proof.
The paper is organized as follows. In Sect. 2, we analyze the local bifurcations, namely pitchfork bifurcation, Hopf bifurcation, Bautin bifurcation and codimension 2 Bogdanov-Takens bifurcation with symmetry. In Sect. 3, we estimate the number of limit cycles in difference parameter regions and curves. In Sect. 4, we study the global bifurcations, namely the different kinds of homoclinic connections and the double limit cycles. In Sect. 5, we give the numerical phase portraits in different parameter regions.

Local bifurcations
Computing the Jacobian matrix at equilibrium E 0 , we have Then, at E 0 the determinant det(J 0 ) = 1 − 1/a and the trace tr(J 0 ) = −bξ , implying that E 0 is a saddle if a < 1, stable focus or node if a > 1 and b > 0, and unstable focus or node if a > 1 and b < 0. When a = 1, equilibrium E 0 is a stable node if b > 0, and unstable node if b < 0 by applying [8,Theorem 2.19]. By the symmetry of vector field (3), equilibrium E L is of the same type as E R . The Jacobian matrix at E R is Calculate det(J R ) = 1 − a 2 and tr(J R ) = −ξ(b + 3 − 3a 2 ). Hence, E R is a stable focus or node if a < 1 and b + 3 − 3a 2 > 0, and unstable focus or node if a < 1 and b + 3 − 3a 2 < 0.

Pitchfork bifurcation
From the expressions of the equilibria E L , E 0 and E R , it follows immediately that a pitchfork bifurcation occurs at the origin of coordinates when a = 1; i.e., for a ≥ 1 we have a unique antisaddle, while for 0 < a < 1 from the previous antisaddle it bifurcates at a = 1 a saddle and two antisaddles. For more details on this kind of bifurcation, see [7,10].

Hopf bifurcations
There are two kinds of Hopf bifurcations, one at the equilibrium E 0 and the other at the equilibria E L and E R , which is essentially the same bifurcation in both, because due to the invariance of system (3) with respect to the symmetry (x, y) → (−x, −y), what occurs at the equilibrium point E L occurs to its symmetric E R .

Hopf bifurcation at E 0
The next result characterizes the Hopf bifurcation at the equilibrium point E 0 , it is proved using the averaging theory, in this way we also can estimated the shape of the limit cycle bifurcating from E 0 , and we avoid the computation of the Liapunov constant.

Proposition 2
The following statements hold for differential system (3).
a. If a > 1, b = 0 and ξ > 0, then a Hopf bifurcation takes place at the equilibrium point located at the origin of coordinates, and the limit cycle γ bifurcating from this equilibrium exists for b < 0 sufficiently small. b. For ε > 0 sufficiently small if b = −βε 2 < 0, then the limit cycle γ passes through the point Moreover, this limit cycle is stable.
Since we want to study the Hopf bifurcation at the origin of coordinates, we blow up the origin doing the scaling r = ε R, differential system (4) taking as new independent variable the θ becomes In order to apply the averaging theory described in appendix, we need that differential equation (5) starts at least with order ε. So we do the change of variables R → ρ defined by ρ.
Then, differential equation (5) in the new variable ρ writes Differential equation (6) is written into normal form (34) for applying the averaging theory summarized in appendix, using the notation of appendix we only need to take n = 1, x = ρ, t = θ, μ = ε 2 , F 1 (t, x) = F 1 (θ, ρ) and T = 2π , and all the necessary assumptions for applying the averaging theory described in appendix are satisfied. Then, we compute Since β > 0 the averaged function f 1 (ρ) has a unique positive zero, ρ = 2 √ β/3 which satisfies the condition In fact, this last expression is positive because a > 1 and ξ > 0, and consequently, by the results described in appendix differential equation (6) has a periodic solution ρ(θ, ε) satisfying that Moreover, this periodic solution ρ(θ, ε) is unstable because derivative (7) is positive. Now we will go back through the changes of variables for obtaining the periodic solution bifurcating from the equilibrium at the origin of coordinates of differential system (3). Thus, the periodic solution ρ(θ, ε) satisfying initial condition (8) in the variables of differential system (5) becomes the periodic solution This periodic solution in differential system (4) becomes (r (t, ε), θ(t, ε)) with and it pass through point in the coordinates (r, θ). Finally, this periodic solution in the coordinates of system (3) is the periodic solution (x(t, ε), y(t, ε)) given by passing through point (9) now in coordinates (x, y). Therefore, when ε → 0 such periodic solution tends to the origin, it is a periodic solution of a Hopf bifurcation. We remark that the periodic solution ρ(θ, ε) was an unstable limit cycle, but due to the fact thatθ is negative in a neighborhood of the origin, when we pass the unstable limit cycle R(θ, ε) to the periodic solution (r (t, ε), θ(t, ε)) it changes to a stable limit cycle.
In short, Proposition 2 shows the existence of the Hopf bifurcation surface H 1 .

Proposition 3
The following statements hold for differential system (3).
We translate the equilibrium point E R to the origin of coordinates doing the change Writing differential system (10) in polar coordinates, we geṫ r = −r cos θ ξr 2 cos 3 θ + 3ξ 1 − a 2 r cos 2 θ − sin θ − r cos θ + 1 − a 2 × 1 − 1 r 2 cos 2 θ + 2 √ 1 − a 2 r cos θ + 1 sin θ − ε 2 ξβr cos 2 θ, θ = −1 + ξ r 2 cos 2 θ + 3 1 − a 2 r cos θ sin θ cos θ Again since we want to study the Hopf bifurcation now at the origin of coordinates, we blow up the origin doing the scaling r = ε R; then, differential system (11) taking as new independent variable the θ writes Again for applying the averaging theory of appendix, we need that differential equation (12) starts at least with order ε. Hence, we do the change of variables R → ρ defined by ρ.
In particular, Proposition 3 shows the existence of the Hopf bifurcation surfaces H 2 and H 3 .

The Bautin bifurcation curve B 0
The standard or classical Hopf bifurcation in a 2dimensional differential system, i.e., that a limit cycle bifurcates from an equilibrium point, takes place in an equilibrium point with purely imaginary eigenvalues which is not a center because the first Liapunov constant at that equilibrium is not zero. These are the Hopf bifurcations studied in Sect. 2.2.2. But when the first Liapunov constant is zero, and also can bifurcate a limit cycle of the equilibrium point if the second Liapunov constant is not zero, such more degenerate Hopf bifurcation is called for some authors a Bautin bifurcation. See for more details about these Hopf bifurcations [11,Chapter 8].

Codimension 2 Bogdanov-Takens bifurcation with symmetry
From [15] a codimension 2 Bogdanov-Takens bifurcation with symmetry happens in system (3) in a neighborhood of the curve a = 1 and b = 0. So, in a neigborhood of the intersection point of P, H 1 , H 2 , H L and DL 1 of the bifurcation diagram in Fig. 1

The dynamics near infinity
In this subsection, we will discuss the qualitative properties of the equilibria at infinity, which describe the behavior of the orbits of system (3) when x and y are sufficiently large. Fig. 3

Proposition 4 As shown in
where dt = z 2 dτ . In this local chart, we only need to study the equilibrium B:(0, 0) of system (18), which corresponds to two equilibria I B + and I B − at infinity of system (3) at the endpoints of the positive and negative y-semiaxes, respectively. By Lemmas 1 and 3 of [13, Chapter 2], we only need to discuss the orbits along characteristic directions of system (18) at B. Applying the polar coordinate changes x = r cos θ and y = r sin θ , system (18) can be written H (θ ) = sin 2 θ cos θ cos 2 θ + a 2 sin 2 θ.
Hence, a necessary condition for θ = θ 0 to be an characteristic direction is G(θ 0 ) = 0, which has exactly two roots 0 and π . Except these two directions, there are no directions along which system (18) has orbits connecting B.
Notice that vector field (18) is symmetric with respect to the v-axis. Thus, we only need to discuss the orbits connecting the origin B of (18) in the half-plane z ≥ 0. We will construct some related open quasisectors to determine how many orbits of (18) connect B in the first and the second quadrants.
Observing that system (18) has four horizontal isoclines: the v-axis, the z-axis and where > 0 is a sufficiently small constant. Set The possible vertical isocline is Obviously, the isocline V is tangent to the v-axis at the origin. Set where σ > 0 is a small constant. Hence, if there exist orbits of system (18) connecting B along the direction of the v-axis in the first and the second quadrants, then near the origin must lie in the sector regions V + BH + or V − BH − if a < 1, and V + BL + or V − BL − if a = 1. The directions of vector field of (18), i.e., the directions of arrows, and the positions of the isoclines are shown in Fig. 4.
Firstly, we consider the case a < 1. We can check Lemma 4 in [14] guarantees that no orbits connect B in V + BV. There are also no orbits connecting B in the interior of is not equal to the slopes of the curves tangent to the v-axis. On the other hand, we compute that (∂/∂v)( 1 (v, z)/ 2 (v, z)) < 0 in the generalized normal sector V BH + of class II, i.e.,ṙ > 0 in V BH + and all positive semiorbits starting from the curves BV and BH + go into V BH + . The definition of generalized normal sectors can be seen in [14,Section 2]. Therefore, there exists a unique orbit leaving from B in V BH + by Lemmas 2 and 5 in [14].
Similarly, in case a = 1, we can also prove that exactly one orbit connects B along the v-axis, which lies in V BL + .

Limit cycles
Lemma 5 Assume that a ≥ 1. System (3) has no limit cycles if b ≥ 0 and a unique limit cycle if b < 0. When b < 0, the following conditions are satisfied.
iv. f (x) and g(x) satisfy the Lipschitz condition in any bounded interval.
Then, by Theorem 4.1 of [17, Chapter 4], system (3) has a unique limit cycle, which is stable. Since the phase portrait of system (3) is symmetric with respect to the point E 0 , the small limit cycles surrounding E L are of the same type as that surrounding E R . Hence, in what follows we only consider the small limit cycles around E R .
Since the proof is similar to Lemma 4 of [5], we omit it. Lemma 6 shows that φ 1 (a, ξ) < a 2 − 1 and the phase portrait (a) in Fig. 2 in Theorem 1 is obtained.
Consider equation where bothF(z) andF (z) are continuous in [0, z 0 ), andF(0) = 0. Let L J denote the integral curve of (19) passing through the point P(z J ,F(z J )) on the curve y =F(z). Also, let y = ϕ J (z) and y =φ J (z) represent the orbit segments of L J below and above the curve y =F(z). When 0 < z < z J , we clearly have ϕ J (z) <F(z) <φ J (z) and ϕ J (z) > 0 >φ J (z). Moreover, we introduce the symbol Then, we have for some a 0 .

Lemma 7 (Lemma 4.5 of [17, Chapter 4]) For equation
Lemma 7 will be applied in the following lemma.

Lemma 8
If 0 < a < 1 and b < a 2 − 1, then system (3) has at most two large limit cycles. Proof Assume that system (3) has at least two large limit cycles surrounding the three equilibria E L , E 0 and E R , and that L 1 and L 2 are the most external limit cycles, where L 2 denotes the outer one. We first consider 3a 2 − 3 ≤ b < a 2 − 1. The corresponding phase portrait is shown in Fig. 5a. By the Bendixson criterium, each L i has two intersection points, denoted by B i and C i (i = 1, 2) with the straight line x = x 0 , where x 0 is the abscissa of the equilibrium E R , as shown in Fig. 5a. By the symmetry of the phase portrait for i = 1, 2. On the arcs A 1 B 1 and A 2 B 2 , let y = y 1 (x) and y = y 2 (x), respectively. In fact, for each i = 1, 2, we have 2 and y i (x) is the function corresponding to the arc A i B i . Thus, Similarly, we obtain that Setting z = x 0 g(s)ds, from (20), we get By Lemma 7, in order to prove the inequality , we only need to prove thatF(z(x 0 )) = 0,F(z) > 0 andF(z)F (z) is nondecreasing for z > z(x 0 ). Clearly, we haveF(z(x 0 )) = 0 andF(z) > 0 for z > z(x 0 ). Note that .
where the three factors of the last line are positive and increasing. Therefore, is positive and increasing. By Lemma 7, we have that Now we consider b < 3a 2 − 3. The corresponding phase portrait is shown in Fig. 5b. By the Bendixson criterium again, each L i has two intersection points with the straight line x = x 0 , being x 0 the unique positive zero of F (x) when x > 0, and the intersection points are denoted by B i and C i (i = 1, 2). Then, inequality (21) can be proved in a similar way to the case Similarly, we obtain that Now, using again Lemma 7, we only need to prove that where all of the factors of the last two lines are positive and increasing. Therefore, is positive and increasing. By Lemma 7, we have and therefore (21) follows. However, it is impossible to have two attracting (repelling) limit cycles surrounding the same equilibrium (equilibria) adjacent one to the other. So, from inequality (21) and the repelling of the infinity, we obtain that system (3) has at most three large limit cycles, where the outer one is stable, the middle one is semistable, the inner one is stable. Clearly, for fixed values a and ξ , system (3) is a family of generalized rotated vector fields with respect to the parameter b. Assume that system (3) has exactly three large limit cycles. By Theorem 3.5 of [17, Chapter 4], the outer limit cycle and the inner one neither split, nor disappear as b varies monotonically. By Theorem 3.4 of [17, Chapter 4], the middle limit cycle will bifurcate into at least one stable and one unstable cycle when b varies in the suitable direction. This is a contradiction. Therefore, system (3) has at most two large limit cycles. If the two large limit cycles exist, we can obtain that the outer limit cycle is stable and the inner one is unstable.

Proposition 9
Consider system (22), which satisfies has at most one solution increasing) for α < x < 0; Then, system (22) has at most one closed orbit in the region {(x, y) ∈ R 2 : α < x < β}. The closed orbit is simple and unstable (resp. stable) if it exists. Proof Assume that system (22) exists a limit cycle γ , as shown in Fig. 6a. In the following, we will ascertain the sign of Clearly, w =F(x) has two inverse functions, x 1 (w) (respectively, x 2 (w)), on the right (resp. left) side of the origin. The functionsλ(x i (w)) will be denoted simply by λ i (w). By w =F(x), we rewrite system (22) intȯ Let y 1 (w) and y 2 (w)(resp. z 1 (w) and z 2 (w)) be functions determined by the orbits of (24) below (resp. above) the line y = w, which correspond the parts of the trajectories of system (22) below (resp. above) the curve y =F(x) and depend whether they are to the left or right of the origin. From condition (iv), λ 1 (w) = λ 2 (w) has at most one root. When the equation λ 1 (w) = λ 2 (w) has no roots, by the comparison theorem we obtain either z 1 (w) > z 2 (w) and y 1 (w) < y 2 (w) or z 1 (w) < z 2 (w) and y 1 (w) > y 2 (w). By Green formula, we have which contradicts γḡ (x)dx + (y −F(x))dy = 0. When the equation λ 1 (w) = λ 2 (w) has a unique root and the equation z 1 (w) = z 2 (w) has no roots, we can similarly prove that (22) has no limit cycles. When the equation λ 1 (w) = λ 2 (w) has a unique root, applying again the comparison theorem the equation z 1 (w) = z 2 (w) has at most one root. Now we only need to discuss the system λ 1 (w) = λ 2 (w) and z 1 (w) = z 2 (w) having a unique root. Here we consider that condition (v) holds. By (w, y) → (μw, μy) with μ := w * 2 /w * 1 > 1, where w * i is the intersection of z i (w) and the line y = w, system (25) deduces First, we consider that λ 1 (w) > λ 2 (w) as w → 0. Moreover, the limit cycle γ of (22) with x > 0 corresponds to the broken curve in Fig. 6b and y 2 (w), z 2 (w) are changed into two new functions, denoted y 2 (w) andz 2 (w), respectively. It is easy to check that y 2 (w) = μỹ 2 (w/μ) and z 2 (w) = μz 2 (w/μ). By μ > 1 and the increase ofλ(w)/w, we getλ(μw)/μ > λ(w). Furthermore, using again the comparison theorem to (25) and (26), we obtain y 1 (w) <ỹ 2 (w) and z 1 (w) >z 2 (w). Thus, So in the region {(x, y) ∈ R 2 : α < x < β} γ is unstable and simple if it exists. Moreover, it is impossible to have two attracting (repelling) limit cycles surrounding the same equilibrium adjacent one to the other. Therefore, the uniqueness has also been proved. For the case λ 1 (w) < λ 2 (w) as w is small, we can prove that γ f (x)dt > 0 in a similar way to the case λ 1 (w) > λ 2 (w). So, in the region {(x, y) ∈ R 2 : α < x < β}, the limit cycle γ is stable and simple if it exists. Therefore, we have completed this proof.
The proof of Proposition 9 gives the following corollary directly. Under the preparations of Proposition 9 and Corollary 10, we obtain the existence of small limit cycles on the parameter surfaces H 2 and H 3 as follows.
Proof If 1/ √ 3 ≤ a < 1 and b = 3a 2 − 3, i.e., on the parameter curve H 2 , we obtain Fig. 7 by Proposition 4 and Lemma 11, which shows the existence of a Poincaré-Bendixson annulus, i.e., any trajectory starting at a point of the boundary curves of the annulus enters (or leaves) the annulus, and inside the annulus there is no equilibrium points. So, the existence of some large limit cycles is obtained. Assume that system (3) has two large limit cycles. Let the outer limit cycle and inner one denoted by γ 2 and γ 1 . Therefore, γ i div(y − F(x), −g(x))dt ≤ 0 for i = 1, 2. By (21), γ 1 div(y − F(x), −g(x))dt > γ 2 div(y − F(x), −g(x))dt. However, it is impossible to have two attracting (repelling) limit cycles surrounding the same equilibrium (equilibria) adjacent one to the other. Therefore, So γ 1 is a semistable limit cycle. Since the vector field of (3) is rotating with respect to b, by Theorem 3.4 of [17, p. 211] there is a stable limit cycleγ 2 near γ 2 for a perturbation of b, and two limit cyclesγ 1 ,γ 1 (γ 1 is smaller thanγ 1 ) near γ 1 . By (21) and stabilities of equilibria, we obtain )dt. Therefore, system (3) has at most one large limit cycle. Thus, the uniqueness of limit cycle is proved and the phase portrait (d) in Fig. 2 in Theorem 1 is obtained.
Moreover, we can get the phase portraits (i), (e) and (n) in Fig. 2 in Theorem 1 by Proposition 3, Lemmas 8, 11, 12 and the continuity of the vector fields.

Global bifurcation
The aim of this section is to show that the global bifurcation surfaces H L, DL 1 , DL 2 for homoclinic loops and double limit cycles in Fig. 1 exist and how they are located in the bifurcation diagram of system (3). In the following results, we discuss the existence and location of large limit cycles and figure-eight homoclinic loops of system (3).
Proof We consider equivalent system (31) of system (3). When b = a 2 −1, system (31) has no limit cycles by Lemma 6 and equilibria E R , E L are stable by the beginning part of Sect. 2. Thus, the separatrices of the saddle at the origin of (31) are shown in Fig. 8a when b = a 2 − 1. When b = 3a 2 − 3 and √ 3/3 < a < 1, system (31) has exactly one large limit cycle; by Lemma 12 and equilibria, E R , E L are unstable by Sect. 2.2.3, yielding that the separatrices of the saddle at the origin of (31) are shown in Fig. 8b. Moreover, (a, b, ξ) = (1, 0, ξ  *  ) is a Bogdanov-Takens bifurcation point of codimension 2 by [15] and (a, b, ξ) = (0, b * , ξ * ) is a grazing bifurcation point by [5], where ξ * > 0 is an arbitrary fixed value and b * < −3 is a certain value. As it was proved in [5], the separatrices of the saddle move monotonically when a, ξ are arbitrarily fixed and b increases. Hence, when a and ξ are fixed, there exists a unique value b = ϕ(a, ξ) such that system (31) has exactly one figure-eight loop. Note that div(y, −g(x) − f (x)y)| (x,y)=(0,0) = −bξ > 0 for b < 0. By [7,Chapter 3], the figure-eight loop is unstable. Moreover, from Fig. 3 and Poincaré-Bendixson annulus Theorem, system (31) has at least one large limit cycle if b = ϕ(a, ξ). Using an analogous argument of Lemma 12, we can show that system (31) has at most one large limit cycle in this case. Therefore, statement (iv) is proved.
For arbitrarily fixed values of a and ξ , system (31) is a rotational family of vector fields (see [17] for definitions and properties) with respect to the parameter b, implying that unstable limit cycles increase and stable ones decrease in size when b increases. Furthermore, the double limit cycle which is stable in its outside part splits into a pair of limit cycles. Hence, when a and ξ are arbitrarily fixed, there exists a unique value b = φ 1 (a, ξ) satisfying that system (31) has exactly one large limit cycle, which is semistable. Thus, statements (ii) and (iii) are proved. By properties of rotated system (31) in b, statement (i) can be proved by a similar way as case (iii).

Lemma 15 The surface DL 2 lies in the region
Proof By Lemma 11, system (3) has no small limit cycles if 1/ √ 3 ≤ a < 1 and b = 3a 2 − 3, and at most one small limit cycle around equilibria E R or E L if 0 < a < 1/ √ 3 and b = 3a 2 − 3. Therefore, DL 2 cannot intersect with the curve b = 3a 2 − 3 except when a = 1/ √ 3. Computing the trace at E 0 ,  {(a, b, ξ) ∈ G : 0 < a < 1, 3a 2 − 3 < b < a 2 − 1}. For fixed a and ξ , we obtain that system (3) has a small semistable limit cycle 0 when b = b 0 from the rotational properties of system (3) with respect to b. Now, given b, ξ and a perturbation a → a + ε, there exists a solutionφ(t, x 0 , y 0 ) for t ∈ (0, T ) which lies in the small neighborhood of 0 from the continuous depen- Fig. 11 Simulations with three equilibria dence of solutions on the parameters and initial conditions, where T is the periodic of 0 and ε > 0 is small. Then, using again the properties of the rotational vector fields, we can take a suitable parameter b 0 + ε 1 such that system (3) has a new semistable small limit cycleˆ 0 , because the homoclinic loops cannot be semistable, when ε 1 > 0 is small. By continuity, system (3) has a small semistable limit cycle in the parameter curve either H 2 , or H 3 , which leads to a contradiction. Therefore, system (3) has at most one small limit cycle in G 4 . Moreover, by Propositions 13 and 14, it follows that DL 2 lies in G 3 . We have completed the proof.
By a similar discussion as that in the proof of Lemma 15, the phase portraits (f)-(h) and (j)-(m) of system (3) in Fig. 2 of Theorem 1 are obtained from the properties of rotational vector fields, continuity, Lemma 15 and the results of Sect. 3.

Remark 1
The existence of surface DL 2 of double small limit cycles can be guaranteed by the Bautin bifurcation, as seen in Sect. 2.2.3. If the graph of DL 2 can be expressed as a function b = φ 2 (a, ξ) and has not any knot, there are at most two small limit cycles around the equilibria E R or E L for system (3), as shown in Fig. 2. If the graph of DL 2 has knots, system (3) will have more than two small limit cycles in the parameter region V in Fig. 1. For example, system (3) has exactly two small limit cycles for special parameter value, where the two ones are quadruple, as shown in Fig. 9.
From [5], we can obtain that a pair of grazing loop is stable for a = 0 if they exist. However, the pair of homoclinic loops have to be unstable if they exist when a = 0. In fact, H L, DL 1 , DL 2 have a common intersection point for the limit value a = 0, i.e., ϕ(0, ξ) = φ 1 (0, ξ) = φ 1 (0, ξ) < −3. Since sys- Fig. 12 Simulations with three equilibria tem (31) is a rotational vector field with respect to the parameter b, the manifolds of E 0 move monotonically as a, ξ are fixed and b increases, see [5,6]. Therefore, it is worthwhile to note that H L, DL 1 and DL 2 have no intersection points except at endpoints. Summarizing the previous results, we can obtain Theorem 1, as shown in Fig. 1.

Numerical examples
In this section, we give several numerical examples of previous results.
Example 1 Let a = √ 2 and ξ = 1. When b = 1, the system has a unique equilibrium (0, 0), which is a sink, and no limit cycles, as shown in Fig. 10a.
However, when b = − 1 the system has a unique equilibrium, the origin (0, 0) which is a source. Fur-thermore, from Lemma 5, there is a unique limit cycle, which is stable, as shown in Fig. 10b.
Example 2 Let ξ = 1. When a = √ 2/2 and b = − 2, the system has three equilibria and exactly one large limit cycle, as shown in Fig. 11a.
When a = 0.3 and b = − 2.77, the system has three equilibria and exactly two small limit cycles, as shown in Fig. 11b.
When a = √ 2/2 and b = −0.5, the system has three equilibria and no limit cycles, as shown in Fig. 11c.
Example 3 Let ξ = 0.1. When a = 0.9 and b = − 0.46, the system has three equilibria and exactly two large limit cycles, as shown in Fig. 12a.
When a = 0.9 and b = − 0.555, the system has three equilibria, exactly two small limit cycles and one large limit cycle, as shown in Fig. 12b. When a = 3 √ 2/10 and b = − 2.461, the system has three equilibria, exactly two small limit cycles and two large limit cycles, as shown in Fig. 12c.
When a = √ 5/5 and b = − 2.41, the system has three equilibria, exactly four small limit cycles and one large limit cycle, as shown in Fig. 12d.

Appendix: Averaging theory of first and second order
The averaging theory of second order for studying specifically periodic orbits can be found in [12] for S 3 differential systems and in [1] for Lipschitz differential systems, see also Chapter 11 of [16]. Here we present a brief summary with the results that we need for studying the Hopf bifurcation of the SD oscillator systems.
Consider the differential systeṁ where F 1 , F 2 : R×D → R, R : R×D×(−μ 0 , μ 0 ) → R are continuous functions, T -periodic in the first variable and D is an open subset of R n . Assume: i. F 1 (t, ·) ∈ S 2 (D), F 2 (t, ·) ∈ S 1 (D) for all t ∈ R, F 1 , F 2 , R, D 2 x F 1 , D x F 2 are locally Lipschitz with respect to x, and R is twice differentiable with respect to μ. We define the functions F k0 : D → R for k = 1, 2 as follows Then, for |μ| > 0 sufficiently small there exists a T -periodic solution x(t, μ) of system (34) such that x(0, μ) → a when μ → 0. If for the i for which f i (a) = 0 the real part of all the eigenvalues of the Jacobian matrix D x ( f i )(a) are negative, then the periodic solution x(t, μ) is asymptotically stable, if some eigenvalue has a positive real part then it is unstable.
The averaging theory of first order takes place when f 1 (x) ≡ 0. If f 1 (x) ≡ 0 and f 2 (x) ≡ 0, we say that that we work with the averaging theory of second order.