This is a preprint of: “Bifurcations from families of periodic solutions in piecewise differential systems”, Jaume Llibre, Douglas D. Novaes, Camila A. B. Rodrigues, Phys. D, vol. 404, 132342:1– 12, 2020. DOI: [10.1016/j.physd.2020.132342] BRANCHING OF LIMIT CYCLES FROM FAMILIES OF PERIODIC SOLUTIONS IN PIECEWISE DIFFERENTIAL SYSTEMS JAUME LLIBRE, DOUGLAS D. NOVAES AND CAMILA A. B. RODRIGUES Abstract.Consider a differential system on the form k X 0 x = F0 (t, x) + εi Fi (t, x) + εk+1 R(t, x, ε), i=1 where Fi : S1×D → Rmand R : S1 × D × (−ε0 , ε0 ) → Rm are piecewise C k+1 functions and T -periodic in the variable t. Assuming that the unperturbed system x0 = F0 (t, x) has a d-dimensional submanifold of periodic solutions with d < m we use the Lyapunov-Schmidt reduction method and the averaging theory to study the existence of limit cycles of the above differential system. Keywords Lyapunov-Schmidt reduction · periodic solution · averaging method · non- smooth differential system · piecewise smooth differential system Mathematics Subject Classification (2000) 34C29 · 34C25 · 37G15 · 34C07 1. Introduction and Statement of the main result 1.1. Introduction. The study of the existence of invariant sets, in especial periodic solutions, is very important for understanding the dynamics of a differential system. A limit cycle of a differential system is a periodic solution isolated in the set of all periodic solutions of the differential system. It is well known that there exists a relation between the periodic solutions of a system and the zeros of a Poincaré map and the Jaume Llibre Departament de Matematiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Catalonia, Spain Email: jllibre@mat.uab.cat Douglas D. Novaes Departamento de Matemática, Universidade Estadual de Campinas, Rua Sérgio Baruque de Holanda, 651, Cidade Universitária Zeferino Vaz, 13083-859, Campinas, São Paulo, Brazil, Email: ddnovaes@ime.unicamp.br Camila A. B. Rodrigues Departamento de Matemática, ICMC-Universidade de São Paulo, 13560-970, São Carlos, São Paulo, Brazil, Email: camilaap@icmc.usp.br 1 http://www.gsd.uab.cat 2 J. LLIBRE, D.D. NOVAES AND C.A.B. RODRIGUES displacement function. In this sense the averaging theory is one of the important tools to detect periodic solutions in m-dimensional systems on the form k X (1) x0 = F0 (t, x) + εi Fi (t, x) + εk+1 R(t, x, ε). i=1 A classical introduction to the averaging theory can be found in [21, 22]. Consider the unperturbed system x0 = F0 (t, x) and its set of initial conditions whose orbits are periodic denoted here by Z. Assume that the set Z is a d-dimensional submanifold of Rm such that either dim(Z) = d = m) or dim(Z) = d < m. When dim(Z) = d = m we have many works that study the number of limit cycles of system (1). Assuming that k ∈ {1, 2}, F0 ≡ 0 and F1 , F2 are T - periodic functions in the first variable and locally Lipschitz in the second variable Buica and Llibre proved in [3] that the number of limit cycles of (1) is controlled by the number of zeros of some functions called average functions that depend on F1 if k = 1 and of F1 and F2 if k = 2. In [7] the authors studied the case where F0 is zero or not and Fi are analytic functions for every k = 1, 2, . . . , n, and in [16] it was studied the averaging theory at any order when the functions Fi are only continuous and T -periodic on the first variable. The averaging theory can be extended to discontinuous differential systems. The study of discontinuous differential systems is important in many fields of the applied sciences because many problems of physics, engineering, economics, and biology are modeled using differential equations with discontinuous right-hand side, see for in- stance [2, 6, 20]. So, there is a natural interest in studying the averaging theory for discontinuous systems. This was the main objective of the works [11, 12, 14, 18]. When dim(Z) = d < m only the averaging theory is not enough to study the number of limit cycles of the systems and other techniques need to be employed together, as the Lyapunov-Schmidt reduction method. In the case that Fi are smooth functions we have the works [4, 5, 8]. If the functions Fi are not smooth or even continuous we have the works [13, 14], where the authors studied some classes of these systems. A piecewise smooth vector field defined on an open bounded set U ⊂ Rm is a function F : U → Rm which is smooth except on a set Σ of zero measure, called the discontinuity set of the vector field F . We suppose that U \ Σ is a finite union of disjoint open sets Ui , i = 1, 2, . . . , n, where the restriction Fi = F U can be extended continuously to Ui . i The orbit of F at a point p ∈ Ui is defined as usual for a differential system. But if p ∈ Σ then the definition of this orbit through p is more delicate. In [9] Filippov used the theory of differential inclusion (see [1]) to give the definition of what is a local orbit at the points of discontinuity where the set Σ is locally a codimension one embedded submanifold of Rm . If p ∈ Σ and Up is a small neighborhood of p then we divide Up \ Σ in two disjoint open sets Up+ and Up− and write F ± (p) = F U ± (p). p In short, let S ⊂ Σ be an embedded hypersurface in S1 × D and Tp S denotes the tangent space of S at the point p. Let l(p) be the segment connecting the vectors F + (p) and F − (p) and the crossing region of the hypersurface S is the set Σc (S) = {p ∈ S : l(p) ∩ Tp S = ∅}. For a point p on the crossing region the local orbit of F at p is given BRANCHING OF LIMIT CYCLES IN PIECEWISE DIFFERENTIAL SYSTEMS 3 as the concatenation of the local trajectories of F ± at p. In this case we say that the orbit crosses the set of discontinuity and that p is a crossing point. When p is not a crossing point we say that p is a sliding point and the local trajectory of F at p slides on Σ. For more details on the Filippov conventions see [9, 10]. In what follows we describe how to use the averaging theory and Lyapunov-Schmidt reduction method for computing isolated periodic solutions of the piecewise smooth differential systems. Then, we set the class of non-autonomous discontinuous piecewise smooth differential equations that we are interested as well as our main result (Theorem A). 1.2. Lyapunov-Schmidt reduction. Consider the function k X (2) g(z, ε) = εi gi (z) + O(εk+1 ), i=0 where gi : D → Rm is a C k+1 function, k ≥ 1, for i = 0, 1, . . . , k, in which D an open bounded subset of Rm . For d < m, let V be an open bounded subset of Rd and β : V → Rm−d a C k+1 function such that (3) Z = {zα = (α, β(α)) : α ∈ V } ⊂ D. The main hypothesis is (Hα ) the function g0 vanishes on the d-dimensional submanifold Z of D. In [5] the authors used the Lyapunov-Schmidt reduction method to develop the bifurcation functions of order i, for i = 0, 1, . . . , k, which for |ε| 6= 0 sufficiently small control the existence of branches of zeros z(ε) of system (2) that bifurcate from z(0) ∈ Z. In this subsection we present the results developed in that work and those that we shall need later on. First we present some notation. Consider the projections onto the first d coordinates and onto the last m − d coordinates denoted by π : Rd × Rm−d → Rd and π ⊥ : Rd × Rm−d → Rm−d , respectively. Also, for a point z ∈ Z we write z = (a, b) ∈ Rd × Rm−d . Let L be a positive integer, let x = (x1 , x2 , . . . , xm ) ∈ D, t ∈ R and yj = (yj1 , . . . , yjm ) ∈ Rm for j = 1, . . . , L. Given G : R × D → Rm a sufficiently smooth function, for each (t, x) ∈ R × D we denote by ∂ L G(t, x) a symmetric L–multilinear map which is applied J to a “product” of L vectors of Rm , which we denote as L j=1 yj ∈ R mL . The definition of this L–multilinear map is L K n X L ∂ L G(t, x) ∂ G(t, x) yj = y1i . . . yLiL . ∂xi1 . . . ∂xiL 1 j=1 i1 ,...,iL =1 We define ∂ 0 as the identity functional. 4 J. LLIBRE, D.D. NOVAES AND C.A.B. RODRIGUES The bifurcation functions fi : V → Rd of order i are defined for i = 0, 1, . . . , k as i X X l K 1 (4) fi (α) = πgi (zα ) + ∂bL πgi−l (zα ) γj (α)cj , c1 ! c2 !2!c2 . . . cl !l!cl l=1 Sl j=1 where the γi : V → Rm−d , for i = 1, 2, . . . , k, are defined recursively as γ1 (α) = −∆−1 ⊥ α π g1 (zα ) and X i−1 K −1 1 I0 ⊥ γi (α) = −i!∆α c2 . . . c ci−1 b ∂ π g0 (z α ) γj (α)cj (5) c1 ! c2 !2! i−1 !(i − 1)! Si0 j=1 i−1 XX K l  1 L ⊥ + c c ∂b π gi−l (z α ) γj (α)cj . c1 ! c2 !2! 2 . . . cl !l! l l=1 Sl j=1 We denote by Sl the set of all l-tuples of non-negative integers (c1 , c2 , . . . , cl ) such that c1 + 2c2 + . . . + lcl = l, L = c1 + c2 + . . . + cl , and by Si0 the set of all (i − 1)-tuples of non-negative integers (c1 , c2 , . . . , ci−1 ) such that c1 + 2c2 + . . . + (i − 1)ci−1 = i, ∂π ⊥ g0 I 0 = c1 + c2 + . . . + ci−1 and ∆α = (zα ). ∂b About the zeros of the function (2) the authors proved in [5] the following result. Theorem 1. Let ∆α denote the lower right corner (m − d) × (m − d) matrix of the Jacobian matrix Dg0 (zα ). Additionally to hypothesis (Hα ) we assume that (i) for each α ∈ V , det ∆α 6= 0; and (ii) f1 = f2 = . . . = fk−1 = 0 and fk is not identically zero. If there exists α∗ ∈ V such that fk (α∗ ) = 0 and det(Dfk (α∗ )) 6= 0, then there exists a branch of zeros z(ε) with g(z(ε), ε) = 0 and |z(ε) − zα∗ | = O(ε). 1.3. The averaging theory. In [5], using Theorem 1, the authors studied high order bifurcation of periodic solutions of the following T -periodic C k+1 with k ≥ 1 differential system k X (6) x0 = F (t, x, ε) = F0 (t, x) + εi Fi (t, x) + O(εk+1 ), (t, z) ∈ S1 × D, i=1 where the prime denotes the derivative with respect to the independent variable t, usually called the time. In their work they assumed that the manifold Z, defined in (3), is such that all solutions of the unperturbed system x0 = F0 (t, x), starting at points of Z are T -periodic and dim Z ≤ m. Consider the variational equation ∂F0 (7) y0 = (t, x(t, z, 0))y, ∂x BRANCHING OF LIMIT CYCLES IN PIECEWISE DIFFERENTIAL SYSTEMS 5 where x(t, z, 0) denotes the solution of system (6) when ε = 0, and we denote a funda- mental matrix of system (7) by Y (t, z). The average function of order i of system (6) is defined as yi (T, z) (8) gi (z) = Y −1 (T, z) , i! where Z t y1 (t, z) = Y (t, z) Y (s, z)−1 F1 (s, x(s, z, 0)) ds, 0 Z t  yi (t, z) = i!Y (t, z) Y (s, z)−1 Fi (s, x(s, z, 0)) 0 (9) X 1 K i−1 I0 + b2 . . . b bi−1 ∂ F0 (s, x(s, z, 0)) yj (s, z)bj b1 ! b2 !2! i−1 !(i − 1)! Si0 j=1 i−1 XX K l  1 L bj + ∂ Fi−l (s, x(s, z, 0)) yj (s, z) ds. b1 ! b2 !2!b2 . . . bl !l!bl l=1 Sl j=1 Using the functions gi stated in (8) are defined the functions fi and γi given by (4) and (5), respectively. Under some assumptions and with Theorem 1 it was proved that the simple zeros of the functions fi provide the existence of isolated periodic solutions of the differential system (6). By a simple zero of a function f we mean a point a such that f (a) = 0 and det(Df (a)) 6= 0, where Df (a) denotes the Jacobian matrix of f at the point a. Remark 1. The functions yi (t, z) could be defined recurrently by an integral equation as done in other works (see [11, 16, 17]). Indeed, we define Z t  y1 (t, z) = F1 (s, x(s, z, 0)) + ∂F0 (s, x(s, z, 0))y1 (s, z) ds, 0 Z t i X X 1 (10) yi (t, z) = i! Fi (s, x(s, z, 0)) + 0 b1 ! b2 !2!b2 . . . bl !l!bl l=1 Sl l K  ·∂ L Fi−l (s, x(s, z, 0)) yj (s, z)bj ds, for i = 2, . . . , k, j=1 and it is not difficult to see that solving this integral equations we obtain the formulae (9). For more details on the results of this subsection 1.2 see [5]. 1.4. Standard form and main result. Let n > 1 be a positive integer. For i = 0, 1, . . . , k and j = 1, 2, . . . , n let Fij : S1 × D → R and Rj : S1 × D × (−ε0 , ε0 ) → R be 6 J. LLIBRE, D.D. NOVAES AND C.A.B. RODRIGUES functions C k+1 where D is an open subset of Rm and S1 ≡ R/(T Z). We define n X Fi (t, x) = χ[tj−1 ,tj ] (t)Fij (t, x), i = 0, 1, ..., k, and j=1 Xn R(t, x, ε) = χ[tj−1 ,tj ] (t)Rj (t, x, ε), j=1 where χA (t) is the characteristic function of A defined as ( 1 if t ∈ A, χA (t) = 0 if t ∈ 6 A. The notation t ∈ S1 ≡ R/(T Z) means that all the above functions are T -periodic in the variable t. Consider the discontinuous and T -periodic differential system k X 0 (11) x = F (t, x, ε) = εi Fi (t, x) + εk+1 R(t, x, ε), i=0 and the submanifold Z of periodic solutions of the unperturbed system (12) x0 = F0 (t, x). The set Σ of discontinuity of system (11) is given by Σ = ({t = 0 ≡ T } ∪ {t = t1 } ∪ . . . ∪ {t = tn−1 }) ∩ S1 × D, where 0 < t1 < t2 < . . . < tn−1 < T . For each j = 1, 2, . . . , n and t ∈ [tj−1 , tj ] we have the differential system k X (13) x0 = εi Fij (t, x) + εk+1 Rj (t, x, ε). i=0 To continue we need to give some definition about system (11). For each z ∈ D and ε sufficiently small we denote by x(·, z, ε) : [0, t(z,ε) ) → Rm the solution of system (11) such that x(0, z, ε) = z, where [0, t(z,ε) ) is the interval of definition for the solution x(t, z, ε). Consider the submanifold Z = {zα = (α, β0 (α)) : α ∈ V̄ }, where V is an open bounded subset of Rm , and β0 : V → Rd−m is a C k function with k ≥ 1. Notice that for each zα ∈ Z, (ti , x(ti , zα , 0) ∈ Σc ), for i ∈ {0, 1, . . . , k}. Indeed, for each j = 1, 2, . . . , n, the set of discontinuity can be locally described by h−1 j (0), where f : S × D → R 1 is hj (t, x) = t − tj . It is known that to show that we are in the crossing region it is sufficient to prove that h∇hj (t, x), F j (t, x)ih∇hj (t, x), F j+1 (t, x)i > 0 (see [16]), where ∇hj (t, x) denotes the gradient vector of the function hj (t, x). Here, ∇hj (t, x) = (1, 0) and h∇hj (t, x), F j (t, x)ih∇hj (t, x), F j+1 (t, x)i = 1 > 0. In [15], the averaging theory was developed assuming dim(Z) = m. Here, we are interested in the case dim(Z) < m. Accordingly, we shall extend the average functions BRANCHING OF LIMIT CYCLES IN PIECEWISE DIFFERENTIAL SYSTEMS 7 (8) and the bifurcation functions (4) obtained in [5] to this class of discontinuous dif- ferential system, providing then sufficient conditions in order to control which periodic solutions of Z, with dim Z = d < m, persists to ε 6= 0 sufficiently small. For system (12) we consider the fundamental matrix Y (t, z) of the variational system ∂ (14) y0 = F0 (t, x(t, z, 0))y, ∂x where Y is an m × m matrix. Notice that, for each j = 1, 2, . . . , n, if xj (t, z, ε) denotes the solution of (13) for tj−1 ≤ t ≤ tj , the function t 7→ (∂xj /∂z)(t, z, 0) is a solution of (14) for tj−1 ≤ t ≤ tj . Recall that the right product of a solution of the variational equation (14) by constant matrix is still a solution of (14). Therefore, the solution Y (t, z) can be built as follows:    Y1 (t, z) if 0 = t0 ≤ t ≤ t1 ,   Y2 (t, z) if t1 ≤ t ≤ t2 , Y (t, z) = .   ..    Yn (t, z) if tn−1 ≤ t ≤ tn = T , where ∂x1 Y1 (t, z) = (t, z, 0), and (15) ∂z  −1 ∂xj ∂xj Yj (t, z) = (t, z, 0) (tj−1 , z, 0) Yj−1 (tj−1 , z), for j = 2, 3, . . . , n. ∂z ∂z Our main result says that the simple zeros of the bifurcation functions (4) also controls the branching of isolated periodic solutions of the nonsmooth system (11). The derivatives ∂ j Fi (s, z), which appears in (4), are computed as follows: n ∂ j Fi X ∂ j Fij (s, z) = χ[tj−1 ,tj ] (s) (s, z). ∂z ∂z j j=1 Theorem A. Let ∆α denote the lower right corner (m − d) × (m − d) matrix of the matrix Id − Y −1 (T, z). We assume that the functions defined by (4) and (8) satisfy f1 = f2 = . . . = fk−1 = 0 and that for each α ∈ V , det(∆α ) 6= 0. If there exists α∗ ∈ V such that fk (α∗ ) = 0, and that det(Dfk (α∗ )) 6= 0, then there exists a T -periodic solution ϕ(t, ε) of (11) such that |ϕ(0, ε) − zα∗ | = O(ε). The paper is organized as follows. In section 2 we present the explicit formulae (8) for the average functions of the nonsmooth differential system (11). In section 3 we prove Theorem A and in Section 4 we give two applications of TheoremA. 2. An algorithm for the bifurcation functions In this section we will provide an algorithm for computing the average functions, defined in (10), for the nonsmooth case. Their expressions are defined recurrently and 8 J. LLIBRE, D.D. NOVAES AND C.A.B. RODRIGUES using Bell polynomials, which can be implemented more easily. This is because soft- wares such as Mathematica and Maple have already implemented these polynomials. In [19], it was proved that the average functions defined for smooth cases can be com- puted using Bell polynomials and in [15] the authors did the same for the nonsmooth case. For each pair of nonnegative integers (p, q), the partial Bell polynomial is defined as X p! Y  xj bj p−q+1 Bp,q (x1 , x2 , . . . , xp−q+1 ) = , b1 !b2 ! . . . bp−q+1 ! j! ep,q S j=1 where Sep,q is the set of all (p − q + 1)-tuple of nonnegative integers (b1 , b2 , . . . , bp−q+1 ) satisfying b1 + 2b2 + . . . + (p − q + 1)bp−q+1 = p, and b1 + b2 + . . . + bp−q+1 = q. Moreover, if g and h are sufficiently smooth functions, then using Bell polynomials we have that l X dl g(h(x)) = g (m) (h(x))Bl,m (h0 (x), h00 (x), . . . , h(l−m+1) (x)). dxl m=1 2.1. Average Functions. In this section we develop a recurrence to compute the average function (8) in the particular case of the discontinuous differential equation (11). So, consider the functions wij : (tj−1 , tj ] × D → R defined recurrently for i = 1, 2, . . . , k and j = 1, 2, . . . , n, as (16) Z t  w11 (t, z) = F11 (s, x(s, z, 0)) + ∂F01 (s, x(s, z, 0))w11 (s, z) ds, 0 Z t wi1 (t, z) = i! Fi1 (s, x(s, z, 0))+ 0 i X X l K  1 · ∂ L Fi−l 1 (s, x(s, z, 0)) 1 wm (s, z)bm ds, b1 ! b2 !2!b2 . . . bl !l!bl l=1 Sl m=1 Z t  wij (t, z) = wij−1 (tj−1 , z) + i! Fij (s, x(s, z, 0))+ tj−1 i X X l K  1 L j j ·∂ Fi−l (s, x(s, z, 0)) wm (s, z)bm ds. b1 ! b2 !2!b2 . . . bl !l!bl l=1 Sl m=1 Since F0 6= 0 the recurrence defined in (16) is an integral equation and the next lemma solves it using Bell polynomials. BRANCHING OF LIMIT CYCLES IN PIECEWISE DIFFERENTIAL SYSTEMS 9 Lemma 2. For i = 1, 2, . . . , k and j = 1, 2, . . . , n the recurrence (16) can be written as follows Z t 1 w1 (t, z) = Y1 (t, z) Y1−1 (s, z)F11 (s, x(s, z, 0))ds, 0  Z t  w1j (t, z) = Yj (t, z) Yj−1 (tj−1 , z)w1j−1 (tj−1 , z) + Yj−1 (s, z)F1j (s, x(s, z, 0))ds , tj−1 Z t  wi1 (t, z) = Y1 (t, z) Y1−1 (s, z) i!Fi1 (s, x(s, z, 0)) 0 i X + ∂ m F01 (s, x(s, z, 0)).Bi,m (w11 , . . . , wi−m+1 1 ), m=2 i−1 X X l i! m 1  + ∂ Fi−l (s, x(s, z, 0)).Bl,m (w11 , . . . , wl−m+1 1 ) ds, l! l=1 m=1 h Z t  wij (t, z) −1 j−1 = Yj (t, z) Yj (tj−1 , z)wi (tj−1 , z) + Yj−1 (s, z) i!Fij (s, x(s, z, 0)) tj−1 i X + j ∂ m F0j (s, x(s, z, 0)).Bi,m (w1j , . . . , wi−m+1 ), m=2 i−1 X X l i! m j  i + j ∂ Fi−l (s, x(s, z, 0)).Bl,m (w1j , . . . , wl−m+1 ) ds. l! l=1 m=1 Proof. The idea of the proof is to relate the integral equations (16) to the Cauchy problem and then solve it. For example, if i = j = 1 the integral equation is equivalent to the following Cauchy problem ∂w11 (t, z) = F11 (t, x(t, z, 0)) + ∂F01 (t, x(t, z, 0)) w11 with w11 (0, z) = 0, ∂t and solving this linear differential equation we get the expression of w11 (t, z) described in the statement of the lemma. For more details see [15].  Now, we provide a formula for the average functions (8) for the class of discontinuous differential systems studied in this paper. Proposition 3. For i = 1, 2, . . . , k, the average function (8) of order i is win (T, z) gi (z) = Yn−1 (T, z) . i! Proof. For each i = 1, 2, . . . , k we define Xn wi (t, z) = χ[tj−1 ,tj ] (t)wij (t, z). j=1 10 J. LLIBRE, D.D. NOVAES AND C.A.B. RODRIGUES Given t ∈ [0, T ] there exists a positive integer k̄ such that t ∈ (tk̄−1 , tk̄ ] and, therefore, wi (t, z) = wik̄ (t, z). By the proof of Proposition 2 of [15] we obtain Z t  w1 (t, z) = F1 (s, x(s, z, 0)) + ∂F0 (s, x(s, z, 0))w1 (s, z) ds, 0 (17) l  K wi (t, z) ·∂ L Fi−l (s, x(s, z, 0)) wm (s, z)bm ds. m=1 Since by Remark 1 we can consider the functions (9) given implicitly, we compute the derivatives in the variable t of the functions (17) and (10) for i = 1, and we see that the functions w1 (t, z) and y1 (t, z) satisfy the same differential equation. Moreover, for each i = 2, . . . , k, the integral equations (10) and (17), which provide respectively yi and wi , are defined by the same recurrence. Then the functions yi and wi satisfy the same differential equations for i = 1, 2, . . . , k, and their initial conditions coincide. Indeed, let i ∈ {1, 2, . . . , k}, since yi (0, z) = 0 and, by (17), wi (0, z) = 0, it follows that the initial conditions are the same. Applying the Existence and Uniqueness Theorem on the solutions of the differential system we get yi (t, z) = wi (t, z), for all i ∈ {1, 2, . . . , k}.  2.2. Bifurcation Functions. In this section we shall write the bifurcation functions (4) and the functions γi (α) given by (5) in terms of Bell polynomials. Claim 1. The bifurcation function (4) is given by i X X l 1 m fi (α) = πgi (zα ) + ∂ πgi−l (zα )Bl,m (γ1 (α), . . . , γl−m+1 (α)), l! b l=1 m=1 where γ1 (α) = −∆−1 ⊥ α π g1 (zα ) and Xi−1 l −1 i! X m ⊥ γi (α) = −∆α ∂b π gi−l (zα )Bl,m (γ1 (α), . . . , γl−m+1 (α)) l! l=0 m=1 i X  m ⊥ + ∂b π g0 (zα )Bi,m (γ1 (α), . . . , γi−m+1 (α)) . m=2 Proof. The expressions (4) was obtained in [5] using the Faá di Bruno’s formula for the L-th derivative of a composite function. This claim follows just by applying the version of the Faá di Bruno’s formula in terms of the Bell polynomials (see [19, 5]).  BRANCHING OF LIMIT CYCLES IN PIECEWISE DIFFERENTIAL SYSTEMS 11 3. Proof of Theorem A For j = 1, 2, . . . , n let ξj (t, t0 , z0 , ε) be the solution of the discontinuous differential system (13) such that ξj (t0 , t0 , z0 , ε) = z0 . Then, we define the recurrence x1 (t, z, ε) = ξ1 (t, 0, z, ε) xj (t, z, ε) = ξj (t, tj−1 , xj−1 (tj−1 , z, ε), ε), j = 2, . . . , n. Since we are working in the cross region it is easy to see that, for |ε| = 6 0 sufficiently small, each xj (t, z, ε) is defined for every t ∈ [tj−1 , tj ]. Therefore x(·, z, ε) : [0, T ] → R is defined as    x1 (t, z, ε) if 0 = t0 ≤ t ≤ t1 ,     x2 (t, z, ε) if t1 ≤ t ≤ t2 ,    ... x(t, z, ε) =   xj (t, z, ε) if tj−1 ≤ t ≤ tj ,     . ..    xn (t, z, ε) if tn−1 ≤ t ≤ tn = T . Notice that x(t, z, ε) is the solution of the differential equation (12) such that x(0, z, ε) = z. Moreover, the following equality hold xj (tj−1 , z, ε) = xj−1 (tj−1 , z, ε), for j = 1, 2, . . . , n. The next lemma expands the solution xj (·, z, ε) around ε = 0. Lemma 4. For j ∈ {1, 2, . . . , n} and tjz > tj , let xj (·, z, ε) : [tj−1 , tj ) be the solution of (13). Then Xk εi j xj (t, z, ε) = xj (t, z, 0) + w (t, z) + O(εk+1 ). i! i i=1 Proof. First, fixed j ∈ {1, 2, . . . , n}, we use the continuity of the solution xj (t, z, ε) and the compactness of the set [tj−1 , tj ] × D × [−ε0 , ε0 ] to get that Z t Rj (s, xj (s, z, ε), ε)ds = O(ε), t ∈ [tj−1 , tj ]. tj−1 Thus, integrating the differential equation (13) from tj−1 to t, we get X k Z t xj (t, z, ε) = xj (tj−1 , z, ε) + εi Fij (s, xj (s, z, ε))ds + O(εk+1 ), and tj−1 (18) Zi=0 t xj (t, z, 0) = xj (tj−1 , z, 0) + F0j (s, xj (s, z, 0))ds. tj−1 By the differentiable dependence of the solutions of a differential system on its pa- rameters the function ε 7→ xj (t, z, ε) is a C k+1 map. Then, the next step is to compute 12 J. LLIBRE, D.D. NOVAES AND C.A.B. RODRIGUES the Taylor expansion of Fij (t, xj (t, z, ε)) around ε = 0 and for this we use the Faá di Bruno’s Formula about the l-th derivative of a composite function, which guarantees that if g and h are sufficiently smooth functions then X l  K b j dl l! (L) (j) g(h(α)) = g (h(α)) h (α) , dαl b1 ! b2 !2!b2 . . . bl !l!bl Sl j=1 where Sl is the set of all l-tuples of non-negative integers (b1 , b2 , . . . , bl ) satisfying b1 + 2b2 + . . . + lbl = l, and L = b1 + b2 + . . . + bl . For each i = 0, 1, ..., k − 1, expanding Fij (s, xj (s, z, ε)) around ε = 0 we get (19) Fij (s, xj (s, z, ε)) = Fij (s, xj (s, z, 0))+ k−i X X l K εl b b ∂ L Fij (s, xj (s, z, 0)) j rm (s, z)bm , b1 ! b2 !2! . . . bl !l! l 2 l=1 Sl m=1 where j ∂m rm (s, z) = xj (s, z, ε) , ∂εm ε=0 and for i = k (20) Fkj (s, xj (s, z, ε)) = Fkj (s, xj (s, z, 0)) + O(ε). Substituting (19) and (20) in (18) we get Z t k X xj (t, z, ε) = xj (tj−1 , z, ε) + εi Fij (s, xj (s, z, 0))ds tj−1 i=0 k−1 X X k−i X 1 + εl+i b1 ! b2 !2!b2 . . . bl !l!bl i=0 l=1 Sl l ! K L j j bm ·∂ Fi (s, xj (s, z, 0)) rm (s, z) ds + O(εk+1 ). m=1 Then, the proof of the lemma ends using the next two claims. Claim 2. For j = 1, 2, . . . , n we have k X εi xj (t, z, ε) = xj (t, z, 0) + rij (t, z) + O(εk+1 ). i! i=1 Claim 3. The equality rij = wij holds for i = 1, 2, . . . , k and j = 1, 2, . . . , n.  Proof of Theorem A. Consider the displacement function (21) h(z, ε) = x(T, z, ε) − z = xn (T, z, ε) − z BRANCHING OF LIMIT CYCLES IN PIECEWISE DIFFERENTIAL SYSTEMS 13 It is easy to see that x(·, z, ε) is a T -periodic solution if and only if h(z, ε) = 0. Moreover, to study the zeros of (21) is equivalent to study the zeros of (22) g(z, ε) = Yn−1 (T, z)h(z, ε). From Lemma 4 we have that k X εi (23) xn (T, z, ε) = xn (T, z, 0) + win (T, z) + O(εk+1 ), i! i=1 for all (t, z) ∈ S1 × D. Replacing (23) in (22) it follows that k ! X εi g(z, ε) = Yn−1 (T, z) xn (T, z, 0) − z + win (T, z) + O(ε k+1 ) i! i=1 k X (24) = Yn−1 (T, z)(xn (T, z, 0) − z) + gi (z) + O(εk+1 ) i=1 k X = gi (z) + O(εk+1 ), i=0 where g0 (z) = Yn−1 (T, z)(xn (T, z, 0) − z). From hypothesis (H) the function g0 (z) vanishes on the submanifold Z, therefore hypothesis (Hα ) holds for the function (24). In order to take the derivative of g0 (z) with respect to the variable z we have the next claim. Claim 4. For every j ∈ {1, 2, . . . , n} ∂xj Yj (tj , z) = (tj , z, 0). ∂z The proof will be done by induction on j. For j = 1 the claim is exactly the definition. Suppose that the claim is valid for j = j0 − 1 and we shall prove it for j = j0 . Since xj (tj−1 , z, ε) = xj−1 (tj−1 , z, ε) for all j = 1, 2, . . . , n we have  −1 ∂xj0 ∂xj0 Yj0 (tj0 , z) = (tj0 , z, 0) (tj0 −1 , z, 0) Yj0 −1 (tj0 −1 , z) ∂z ∂z  −1 ∂xj0 ∂xj0 −1 ∂xj0 −1 = (tj0 , z, 0) (tj0 −1 , z, 0) (tj0 −1 , z, 0) ∂z ∂z ∂z ∂xj0 = (tj0 , z, 0). ∂z Hence if z ∈ Z then   ∂g0 −1 ∂x (z) = Y (T, z) (T, z, 0) − Id ∂z ∂z = Y −1 (T, z)(Y (T, z) − Id) = Id − Y −1 (T, z), 14 J. LLIBRE, D.D. NOVAES AND C.A.B. RODRIGUES which has by assumption its lower right corner (m−d)×(m−d) matrix ∆α nonsingular. From here, the result follows from Proposition 3 and Theorem 1.  4. Examples This section is devoted to present some applications of Theorem A. The first one is as 3D piecewise smooth system for which the plane y = 0 is the discontinuous manifold and admits a surface z = f (x, y) foliated by periodic solutions. The second one is a 3D piecewise smooth system for which the algebraic variety xy = 0 is the discontinuous set and the plane z = 0 has a piecewise constant center. For these systems, we compute some of the bifurcations functions in order to study the persistence of periodic solutions. 4.1. Nonsmooth perturbation of a 3D system. Let f : R2 → R and g : R2 → R be differential functions such that g(x, y) = f (x, y) + x∂y f (x, y) − y∂x f (x, y). Consider the nonsmooth vector field     Xε+ (x, y, z), y > 0 (25) Xε (x, y, z) =    Xε− (x, y, z), y < 0 where  Xε+ (x, y, z) = −y + ε(a0 + a1 z) + ε2 (a2 + a3 z), x, −z + g(x, y) , and  Xε− (x, y, z) = −y, x + εb1 z + ε2 (b2 + b3 )z, −z + g(x, y) , with a0 , a1 , a2 , b1 , b2 , b3 ∈ R. Denote the discontinuous se by Σ = {(x, y, z) ∈ R3 : y = 0}. Notice that the surface z = g(x, y) is an invariant set of the unperturbed vector field X0 . Indeed, considering the function fˆ(x, y, z) = z − f (x, y), we get h∇fˆ(x, y, z), X0 (x, y, z)i z=f (x,y) = 0.  Moreover, since X0 (x, y, f (x, y)) = − y, x, x∂y f (x, y) − y∂x f (x, y) we conclude that the invariant set z = f (x, y) is foliate by periodic solutions. Next result gives suficient conditions in order to guarantee the persistence of some periodic solution. Consider the function Z π Z 2π (26) f1 (r) = a1 f (r cos φ, r sin φ) cos φdφ + b1 f (r cos φ, r sin φ) sin φdφ. 0 π Theorem 5. Consider the piecewise vector field (25). Then, for each r∗ > 0, such that f1 (r∗ ) = 0 and f10 (r∗ ) 6= 0, there exists a crossing limit cycle ϕ(t, ε) of X of period Tε = 2π + O(ε) such that ϕ(t, ε) = (x∗ , y ∗ , f (x∗ , y ∗ )) + O(ε) with |(x∗ , z ∗ )| = r∗ . BRANCHING OF LIMIT CYCLES IN PIECEWISE DIFFERENTIAL SYSTEMS 15 In order to apply Theorem A for proving Theorem 5 we need to write system (25) in the standard form. Considering cylindrical coordinates x = r cos θ, y = r sin θ, z = z, the set of discontinuity becomes Σ = {θ = 0} ∪ {θ = t1 } with t0 = 0, t1 = π and t2 = 2π. The differential system (ẋ, ẏ, ż) = Xε+ (x, y, z) in cylindrical coordinates writes r0 (t) = ε(a0 + a1 z) cos θ + ε2 (a2 + a3 z) cos θ, z 0 (t) = g(r cos θ, r sin θ) − z, (a0 + a1 z) sin θ (a2 + a3 z) sin θ θ0 (t) = 1 − ε − ε2 , r r and the differential system (ẋ, ẏ, ż) = Xε− (x, y, z) becomes r0 (t) = εb1 z sin θ + ε2 (b2 + b3 z) sin θ, (27) z 0 (t) = g(r cos θ, r sin θ) − z, b1 z cos θ (a2 + a3 z) cos θ θ0 (t) = 1 + ε + ε2 . r r Notice that, for each j = 1, 2 and tj−1 ≤ θ ≤ tj , we have θ̇(t) 6= 0 for |ε| = 6 0 sufficiently small. Thus, in a sufficiently small neighborhood of the origin we can take θ as the new independent time variable. Accordingly, system (27) becomes r0 (t) ṙ(θ) = = F01 (θ, r, z) + εF11 (θ, r, z) + ε2 F2 (θ, r, z) + O1 (ε3 ), θ0 (t) z 0 (t) ż(θ) = 0 = F02 (θ, r, z) + εF12 (θ, r, z)) + ε2 F22 (θ, r, z) + O2 (ε3 ). θ (t) Considering the notation of Theorem A we have Fi (θ, r, z) = (Fi1 (θ, r, z), Fi2 (θ, r, z)) for each i ∈ {1, 2}. Moreover, for each i ∈ {1, 2} the function Fi (θ, r, z) is written in P the form Fi (θ, r, z) = 2j=1 χ[tj−1 ,tj ] (θ)Fij (θ, r, z). Defining f˜(θ, r) = f (r cos θ, r sin θ) and g̃(θ, r) = g(r cos θ, r sin θ) we write explicitly the expressions of F0 , F1j and F2j for j ∈ {1, 2}, F0 (θ, r, z) = (0, g̃(θ, r) − z),   1 (a0 + a1 z) sin θ F1 (θ, r, z) = (a0 + a1 z) cos θ, (g̃(θ, r) − z) , r   b1 z cos θ F12 (θ, r, z) = b1 z sin θ, − (g̃(θ, r) − z) , r  (a0 + a1 z)2 sin θ cos θ sin θ F21 (θ, r, z) = (a2 + a3 z) cos θ + , 2 (a0 + a1 z)2 sin θ r r  +(a2 + a3 z)r) (g̃(θ, r) − z) ,   2 b21 z 2 sin θ cos θ cos θ 2  F2 (θ, r, z) = (b2 + b3 z) sin θ − , 2 b1 z cos θ − (b2 + b3 z)r (g̃(θ, r) − z) . r r 16 J. LLIBRE, D.D. NOVAES AND C.A.B. RODRIGUES The unperturbed systems is smooth and its solution (r(θ, r0 , z0 ), z(θ, r0 , z0 )) with initial condition (r0 , z0 ) is given by  Z θ  −θ s (28) r(θ) = r(θ, r0 , z0 ) = r0 , z(θ) = z(θ, r0 , z0 ) = e z0 + e g̃(s, r0 )ds . 0 Consequently, a fundamental matrix solution of (14) is given by   ∂(r, z)  1 0  Y (θ, r0 , z0 ) = (θ, r0 , z0 ) =   ,  ∂(r0 , z0 ) G(θ, r0 ) e−θ where G(θ, r0 ) is the derivative of z(θ, r0 , z0 ) with respect to the variable r0 . Notice that, from (28), G(θ, r0 ) does not depend on z0 . Let ε0 > 0 be a real positive number and consider the set Z ⊂ R2 such that Z = {(r, f˜(0, r)) : r > ε0 }. Notice that for (r0 , z0 ) = (r0 , f˜(0, r0 )) ∈ Z we have z(θ, r0 , z0 ) = f˜(θ, r0 ) = f (r0 cos θ, r0 sin θ). Indeed, let w(θ) = f (r0 cos θ, r0 sin θ). So w0 (θ) = ∂x f (r0 cos θ, r0 sin θ)(−r0 sin θ) + ∂y f (r0 cos θ, r0 sin θ)(r0 cos θ) = g(r0 cos θ, r0 sin θ) − f (r0 cos θ, r0 sin θ) = g(r0 cos θ, r0 sin θ) − w(θ) = g̃(θ, r0 ) − w(θ). The second equality holds because g(x, y) = f (x, y) + x∂y f (x, y) − y∂x f (x, y). Hence, for (r0 , z0 ) ∈ Z the solution z(θ, r0 , z0 ) is 2π-periodic. Moreover,    0 0  Id − Y −1 (2π, r, z) =   .  ? 1−e 2π Consequently, ∆α = 1 − e2π 6= 0. Accordingly, all the hypotheses of Theorem A are satisfied. Proof of Theorem 5. Denote by (r, zr ) a point in Z, that is zr = f˜(0, r). Notice that the bifurcation function of first order is f1 (r) = πg1 (r, zr ), where g1 is defined in (8). ∂πg0 Indeed, from definition f1 (r) = πg1 (r, zr ) + (r, zr )γ1 (r). But ∂b g0 (r, z) = Y −1 (2π, r, z)((r, z(2π, r, z))) − (r, z(0, r, z))) = (0, ?), and then πg0 ≡ 0. Moreover, BRANCHING OF LIMIT CYCLES IN PIECEWISE DIFFERENTIAL SYSTEMS 17 Z θ  Z θ  w11 (θ, r, z) = a0 sin θ + a1 z(φ) cos φdφ − z(φ) cos φdφ, G(θ, r) a0 sin θ + a1 0 0 Z θ  ! −θ φ eφ (g̃(φ, r) − z(φ))(a0 + a1 z(φ)) e e G(φ, r)(a0 + a1 z(φ)) cos φ + sin φ dφ , 0 r " Z # θ −1 −1 w12 (θ, r, z) =Y (θ, r, z) Y (π, r, z)w11 (π, r, z) + Y (φ, r, z)F12 (φ, r(φ), z(φ))dφ π Z π Z θ =Y (θ, r, z) a1 z(φ) cos φdφ + b1 z(φ) sin φdφ, 0 π Z π eφ ((a0 + a1 z(φ))(sin φ(g(r cos φ, r sin φ) − z(φ)) − r cos φ G(φ, r)) dφ 0 r Z θ ! b1 eφ z(φ)(cos φ(g(r cos φ, r sin φ) − z(φ)) + r sin φ G(φ, r)) + − dφ . π r Since g1 (r, z) = Y −1 (2π, r, z)w12 (2π, r, z) and f1 (r) = πg1 (r, zr ) it follows that Z π Z 2π (29) f1 (r) = a1 f (r cos φ, r sin φ) cos φdφ + b1 f (r cos φ, r sin φ) sin φdφ. 0 π So, from Theorem A, each positive simple zero of (26) provides an isolated periodic solution of system (25). This concludes this proof.  The next result is an application of Theorem 5. We shall use in its statement the concept of Bessel functions, which are defined as the canonical solutions y(x) of Bessel’s differential equation d2 y dy x2 2 +x + (x2 − α2 )y = 0, α ∈ C. dx dx This equation has two linearly independent solutions. Using Frobenius’ method we obtain one of these solutions, which is called a Bessel function of the first kind, and is denoted by Jα (x). More details about this function can be found in [23]. Corollary 6. Consider the piecewise vector field (25). (a) If f (x, y) = cos x, then the piecewise smooth vector field X admits a sequence of limit cycles ϕi (t, ε) of X of period Tε such that Tε = 2π + O(ε), ϕn (t, ε) = (x∗n , yn∗ , cos(x∗n )) + O(ε), and |(x∗n , zn∗ )| = nπ/2. (b) If f (x, y) = sin x, then the piecewise smooth vector field X admits a sequence of limit cycles ϕi (t, ε) of X of period Tε such that Tε = 2π + O(ε), ϕi (t, ε) = (x∗n , yn∗ , sin(x∗n )) + O(ε), and |(x∗n , zn∗ )| = rn∗ , where each rn is a zero of the Bessel Function of First Kind, J1 (r). 18 J. LLIBRE, D.D. NOVAES AND C.A.B. RODRIGUES Proof. For f (x, y) = cos x, the bifurcation function (29) reads f1 (r) = −(2b1 sin r)/r, and for f (x, y) = cos(x), the bifurcation function (29) reads f1 (r) = a1 πJ1 (r). There- fore the result follows directly from Theorem 5.  Notice that Theorem 5 cannot be applied when f1 is identically zero, which is the case when f (x, y) = 2x2 − y 2 for instance. For these cases we define the function (30) Z π  Z s f2 (r) = a1 cos s G(s, r) cos φ(a0 + a1 f˜(φ, r))dφ 0 0 Z s − e−s eφ (a0 + a1 f˜(φ, r))(r cos φ G(φ, r) + (f˜(φ, r) − g̃(φ, r))dφ 0 ! sin s + a2 + a3 f˜(φ, r) + (a0 + a1 f˜(s, r))2 ds r Z π e−2π (1 + eπ ) + π (a1 e − b1 ) eφ G(φ, r) cos φ(a0 + a1 f˜(φ, r))dφ 2(1 − e2π ) 0 Z π φ e sin φ + (a0 + a1 f˜(φ, r))(g̃(φ, r) − f˜(φ, r))dφ 0 r Z 2π Z  φ ˜ b1 2π φ ˜ + b1 e G(φ, r) sin φf (φ, r)dφ + e cos φ(g̃(φ, r) − f (φ, r))dφ π r π Z 2π 2 + (−b21 cos s(f˜(s, r))2 + sin s(b2 + b3 f˜(s, r))) π r  Z π Z s + 2b1 sin s G(s, r) cos φ(a0 + a1 f˜(φ, r)) + b1 G(s, r) sin φf˜(φ, r)dφ 0 π Z π eφ sin φ +e −s −eφ cos φ G(φ, r)(a0 + a1 f˜(φ, r)) + (g̃(φ, r) − f˜(φ, r))dφ 0 r Z s   ! φ cos φ ˜ + b1 e (f (φ, r) − g̃(φ, r)) − G(φ, r) sin φ dφ ds. π r Theorem 7. Consider the piecewise vector field (25). Assume that f1 ≡ 0. Then, for each r∗ > 0, such that f2 (r∗ ) = 0 and f20 (r∗ ) 6= 0, there exists a crossing limit cycle ϕ(t, ε) of X of period Tε such that Tε = 2π + O(ε), ϕ(t, ε) = (x∗ , y ∗ , f (x∗ , y ∗ )) + O(ε), and |(x∗ , z ∗ )| = r∗ . Proof. As we saw before πg0 ≡ 0. So, from (4), we compute the bifurcation function of order 2 as ∂πg1 (31) f2 (r) = (r, zr )γ1 (r) + πg2 (r, zr ), ∂b BRANCHING OF LIMIT CYCLES IN PIECEWISE DIFFERENTIAL SYSTEMS 19 1 where γ1 (r) = − π ⊥ g1 (r, zr ) and 1 − e2π Z π φ e ((a0 + a1 f˜(φ, r))(sin φ(g(r cos φ, r sin φ) − f˜(φ, r)) − r cos φG(φ, r, z)) π ⊥ g1 (r, zr ) = dφ 0 r Z 2π eφ f˜(φ, r)(cos φ(g(r cos φ, r sin φ) − f˜(φ, r)) + r sin φG(φ, r, z)) − b1 dφ. π r From Proposition 3, we have g2 (r, zr ) = Y −1 (2π, r, z)w22 (2π, r, z)/2, where wij (2π, r, z) is given in Lemma 2. All these functions may be computed to get (31) as (30). Again, from Theorem A, each positive simple zero of (30) provides an isolated periodic solution of system (25). This concludes this proof.  The next result is an application of Theorem 7. Corollary 8. Consider the piecewise vector field (25) and let f (x, y) = 2x2 − y 2 . Assuming a21 + b21 6= 0 define −80b2 (1 − eπ ) A0 = 2 2 1 2 , (1 + eπ )4(15a π 1 b1 − b1 − 14a1 ) − 5π(1 − e )(b1 + 10a1 ) (32) 40a0 ((1 + eπ )(b1 − a1) − a1π(1 − eπ ) A1 = , (1 + eπ )4(15a1 b1 − b21 − 14a21 ) − 5π(1 − eπ )(b11 + 10a21 ) and D = −4A31 − 27A20 . (i) If D > 0 then the piecewise smooth vector field admits at least one limit cycle. Moreover, if A1 < 0 and A0 > 0, then the piecewise smooth vector field admits at least two limit cycles; (ii) If D ≤ 0 and A0 < 0, then the piecewise smooth vector field admits at least one limit cycle. Moreover, in both cases we have a limit cycle ϕ(t, ε) of X of period Tε such that Tε = 2π + O(ε), ϕ(t, ε) = (x∗n , yn∗ , 2(x∗n )2 − (yn∗ )2 ) + O(ε), and |(x∗n , zn∗ )| = rn∗ . Proof. For f (x, y) = 2x2 − y 2 the bifurcation function (30) becomes (33) a0 ((eπ (1 − π) + 1 + π) a1 − (1 + eπ ) b1 ) f2 (r) = −2b2 + r eπ − 1  − (eπ (56 − 50π) + 56 + 50π) a21 + 60 (1 + eπ ) a1 b1 − (eπ (4 − 5π) + 4 + 5π) b21 3 + r . 40 (eπ − 1) . Dividing f2 by a21 + b21 6= 0, we see that the equation f2 (r) = 0 is equivalent to f˜2 (r) = A0 + A1 r + r3 = 0, where A0 and A1 are given in (32). Notice that f˜2 (r) is a polynomial function of degree 3, so it has at least one real root and can be written as f˜2 (r) = r3 − (r1 + r2 + r3 )r2 + (r1 r2 + r1 r3 + r2 r3 )r − r1 r2 r3 , where ri , i = 1, 2, 3 are the zeros of the polynomial. Moreover, the sign of its discriminant D = −4A31 − 27A20 carries information about its number of real roots. 20 J. LLIBRE, D.D. NOVAES AND C.A.B. RODRIGUES If D > 0 the polynomial f˜2 (r) has three simple real roots r1 , r2 and r3 . Since the polynomial has no quadratic term, it follows that r1 + r2 + r3 = 0 and then at least one of these roots must be positive. Moreover, if A1 < 0 and A0 > 0 then there are two changes of sign between the terms of the polynomial and then by Descartes Sign Theorem we get the two positive roots. If D ≤ 0 then there is a pair of complex roots or a double real root. In both cases the condition A0 < 0 implies that at least one root is positive. Now, from Theorem A, each positive simple zero of (33) provides an isolated periodic solution of system (25). This concludes this proof.  4.2. Nonsmooth perturbation of a nonsmooth center. In this example we con- sider a discontinuous differential system in R3 defined in 4 zones (n = 4). Consider the nonsmooth vector field   X1 (u, v, w) if u > 0 and v > 0,  X (u, v, w) if u < 0 2 and v > 0, (34) X(u, v, w) =  X3 (u, v, w) if u < 0 and v < 0,   X4 (u, v, w) if u > 0 and v < 0, where X1 (u, v, w) = (−1 + ε(a1 x + b1 ), 1, −w + ε(c1 x + d1 )) , X2 (u, v, w) = (−1 + ε(a2 x + b2 ), −1, −w + ε(c2 x + d2 )) , X3 (u, v, w) = (1 + ε(a3 x + b3 ), −1, −w + ε(c3 x + d3 )) , X4 (u, v, w) = (1 + ε(a4 x + b4 ), 1, −w + ε(c4 x + d4 )) , with aj , bj , cj , dj ∈ R for all j. Writing in cylindrical coordinates u = r cos θ, v = r sin θ, w = w, the set of disconti- nuity is Σ = {θ = 0} ∪ {θ = t1 } ∪ {θ = t2 } ∪ {θ = t3 } with t0 = 0, t1 = π/2, t2 = π, t3 = 3π/2 and t4 = 2π. For each j = 1, 2, 3, 4 the differential system (u̇, v̇, ẇ) = Xj (u, v, w) in cylindrical coordinates writes k X r0 (t) = gj (θ) + εi (aij r cos2 θ + bij cos θ), i=1 k X w0 (t) = −w + εi (cij r cos θ + dij cos θ), i=1 k ! 1 X θ0 (t) = gbj (θ) − εi (aij r cos θ sin θ + bij sin θ) , r i=1 BRANCHING OF LIMIT CYCLES IN PIECEWISE DIFFERENTIAL SYSTEMS 21 where g1 (θ) = sin θ − cos θ, gb1 (θ) = sin θ + cos θ, g2 (θ) = −(sin θ + cos θ), gb2 (θ) = sin θ − cos θ, g3 (θ) = − sin θ + cos θ, gb3 (θ) = −(sin θ + cos θ), g4 (θ) = sin θ + cos θ, gb4 (θ) = − sin θ + cos θ. Notice that, for each j = 1, 2, 3, 4 and tj−1 ≤ θ ≤ tj , θ̇(t) 6= 0 for |ε| sufficiently small. Thus, in a sufficiently small neighborhood of the origin we can take θ as the new independent time variable by doing r0 (θ) = ṙ(t)/θ̇(t) and w0 (θ) = ẇ(t)/θ̇(t). Taking θ as the new independent time variable we have j j r0 (θ) = F01 (θ, z) + εF11 (θ, z) + O1 (ε2 ), (35) j j w0 (θ) = F02 (θ, z) + εF12 (θ, z) + O2 (ε2 ). Here, z = (r, w) and the prime denotes the derivative with respect to θ. The expressions j j of F01 and F02 for j = 1, 2, 3, 4 are given by 1 r(sin θ − cos θ) 1 −rw 2 r(sin θ + cos θ) 2 rw F01 = , F02 = , F01 = , F02 = , sin θ + cos θ sin θ + cos θ cos θ − sin θ cos θ − sin θ 3 r(sin θ − cos θ) 3 rw 4 r(sin θ + cos θ) 4 −rw F01 = , F02 = , F01 = , F02 = . sin θ + cos θ sin θ + cos θ cos θ − sin θ cos θ − sin θ j j The expressions of F11 and F12 for j = 1, 2, 3, 4 are also easily computed. Nevertheless, we shall omit these expressions because of their size. For each j ∈ {1, 2, 3, 4}, the differential system (35) is 2π-periodic in the variable θ and is written in the standard form with   Fij (θ, z) = Fi1 j j (θ, z), Fi2 (θ, z) , for i = 0, 1. Now, for each j ∈ {1, 2, 3, 4} we compute the solution xj (θ, z, 0) of the unperturbed system j j ṙ(θ) = F01 (θ, z), ẇ(θ) = F02 (θ, z). and this solution is   r − sinrθ+cos sin θ x1 (θ, z, 0) = sin θ+cos θ , we θ ,   −r − cosrθ−sin sin θ −2r x2 (θ, z, 0) = cos θ−sin θ , we θ ,   −r − sinrθ+cos sin θ −2r x3 (θ, z, 0) = sin θ+cos θ , we θ ,   r − cosrθ−sin sin θ −4r x4 (θ, z, 0) = cos θ−sin θ , we θ . 22 J. LLIBRE, D.D. NOVAES AND C.A.B. RODRIGUES We note that in each quadrant the denominators of these four solutions never vanish. Let 0 < r0 < r1 be positive real numbers and consider the set Z ⊂ R2 such that Z = {(α, 0) : r0 < α < r1 }. The solution x(θ, z, 0) of the unperturbed system x0 (θ) = F0 (θ, z) satisfies x(θ, z, 0) = xj (θ, z, 0), for θ ∈ [tj−1 , tj ], and x(2π, z, 0) − x(0, z, 0) = (0, z(1 − e−4r )). Consequently, for each zα ∈ Z, the solution x(θ, z, 0) is 2π-periodic and system (34) satisfies hypothesis (H). Moreover, the fundamental matrix Y (θ, z) is given by    Y1 (θ, z) if 0 = t0 ≤ θ ≤ π/2,  Y (θ, z) if π/2 ≤ θ ≤ π, 2 Y (θ, z) =   Y3 (θ, z) if π ≤ θ ≤ 3π/2,   Y4 (θ, z) if 3π/2 ≤ θ ≤ 2π, where Yj (t, z) are defined by (15). So   1  g4 (θ) 0  Y1 (θ, z) =   − r sin θ ,  e g4 (θ) w sin θ − rg sin θ − g4 (θ) e 4 (θ)   1  g3 (θ) 0  Y4 (θ, z) =   − r sin θ −4r .  e g3 (θ) w(sin θ+4g 3 (θ)) − rg sin θ −4r − g3 (θ) e 3 (θ) Hence,    0 0  Y1 (0, z)−1 − Y4 (2π, z)−1 =   ,  −4w 1 − e4r and then det(∆α ) = 1 − e4r 6= 0 if zα = (α, 0) ∈ Z. Thus, we can compute the bifurcation functions (4) for system (34). For doing this we first obtain the functions (16) corresponding to this system, g0 (θ, z) = (0, w(1 − e4r )), 1 w14 (2π, z) = r(r(a1 + a2 + a3 + a4 ) + 2(b1 − b2 − b3 + b4 )), 2 1 −4r e (−r2 w(6a1 + 3a2 + 2a3 ) − 3r(w(4b1 − 2b2 − b3 ) 3 +e2r (−e2r c4 + c2 + c3 ) + c1 ) + 3(er − 1)(er (c2 + d2 )  +e2r (c3 − d3 ) + e3r (d4 − c4 ) + c1 + d1 )) , and (36) g1 (z) = Y4 (2π, z)−1 w14 (2π, z). BRANCHING OF LIMIT CYCLES IN PIECEWISE DIFFERENTIAL SYSTEMS 23 So, the bifurcation function (4) corresponding to the function (36) becomes 1 f1 (α) = α(α(a1 + a2 + a3 + a4 ) + 2(b1 − b2 − b3 + b4 )), 2 which has a simple zero α∗ . So, from Theorem A, we get the existence of an isolated periodic solution of system (35) for ε sufficiently small. Acknowledgements JL is partially supported by the Ministerio de Economı́a, Industria y Competitividad, Agencia Estatal de Investigación grant MTM2016-77278-P (FEDER), the Agència de Gestió d’Ajuts Universitaris i de Recerca grant 2017 SGR 1617, the European project Dynamics-H2020-MSCA-RISE-2017-777911, and by the European Community grants FP7-PEOPLE-2012-IRSES-316338 and FP7-PEOPLE-2012-IRSES-318999. DDN is partially supported by a FAPESP grant 2016/11471-2 and by the European Commu- nity grants FP7-PEOPLE-2012-IRSES-316338 and FP7-PEOPLE-2012-IRSES-318999. CABR is partially supported by a Ph.D. CNPq Fellowship number 140292/2017-9. References [1] J. P. Aubin and A. Cellina, Differential Inclusions, Springer-Verlag, Berlin, 1984. [2] M. di Bernardo, C. J. Budd, A. R. Champneys, and P. Kowalczyk, Piecewise-Smooth Dynamical Systems: Theory and Applications, Springer, 2008. [3] A. Buica, and J. Llibre, Averaging methods for finding periodic orbits via Brouwer degree, Bull. Sci. Math, 128 (2004), 7–22. [4] A. Buica, J. Françoise and J. Llibre, Periodic solution for non nonlinear periodic differential systems with a small parameter, Comunications on Pure and Applied Analysis 6, (2007), 103–111. [5] M. Cândido, J. Llibre and D.D. Novaes, Persistence of periodic solutions for higher order perturbed differential systems via Lyapunov-Schmidt reduction, Nonlinearity 30 (2017), 3560–3586. [6] S. Coombes, Neuronal networks with gap junctions: A study of piecewise linear planar neuron models, SIAM J. Appl. Math. 7 (2008), 1101–1129. [7] J. Gine, M. Grau, and J. Llibre Averaging theory at any order for computing periodic orbits , Physica D. 250 (2013), 58–65. [8] J. Gine, J. Llibre, K. Wu andX. Zhang Averaging methods of arbitrary order, periodic solutions and integrability, Journal of Differential Equations. 260 (2016), 4130–4156. [9] A. F. Filippov, Differential Equations with Discontinuous Righthand Side, Mathematics and Its Applications, Kluver Academic Publishers, Dordrecht, 1988. [10] M. Guardia, T.M. Seara, M.A. Teixeira, Generic Bifurcations of low codimension of Planar Filippov Systems , J. Differential Equations 250 (2011), 1967–2023. [11] J. Itikawa, J. Llibre and D. D. Novaes, A new result on averaging theory for a classe of discontinuous planar differential systems with applications, Rev. Mat. Iberoam. 33 (2017), 1247– 1265. [12] J. Llibre, A.C. Mereu and D.D. Novaes, Averaging theory for discontinuous piecewise differ- ential systems, J. Differential Equation 258 (2015), 4007–4032. [13] J. Llibre and D. D. Novaes, Improving the averaging theory for computing periodic solutions of the differential equations, Z. Angew. Math. Phys. ZAMP, 66 (2015), 1401–1412. [14] J. Llibre and D.D. Novaes, On the periodic solutions of discontinuous piecewise differential systems, preprint, 2016, arXiv:1504.03008. 24 J. LLIBRE, D.D. NOVAES AND C.A.B. RODRIGUES [15] J. Llibre, D. D. Novaes and C.A.B. Rodrigues, Averaging theory of any order for computing limit cycles of discontinuous piecewise differential systems with many zones, Physica D: Nonlinear Phenomena (2017), 353–354. [16] J. Llibre, D.D. Novaes and M. A. Teixeira, Higher order averaging theory for finding periodic solutions via Brouwer degree, Nonlinearity 27 (2014), 563–583. [17] J. Llibre, D.D. Novaes and M. A. Teixeira, Corrigendum: Higher order averaging theory for finding periodic solutions via Brouwer degree (2014 Nonlinearity 27 563), Nonlinearity 27 (2014), 2417. [18] J. Llibre, D.D. Novaes and M.A. Teixeira, On the birth of limit cycles for non-smooth dy- namical systems, Bull. Sci. Math. 139 (2015), 229–244. [19] D.D. Novaes, An Equivalent Formulation of the Averaged Functions via Bell Polynomi- als., In Extended Abstracts Spring 2016: Nonsmooth Dynamics. Trends in Mathematics 8, Birkhäuser/Springer, Cham, 2017, 141–145 [20] Various, Special issue on dynamics and bifurcations of nonsmooth systems, Phys. D 241 (2012), 1825–2082. [21] J. A. Sanders, F. Verhulst and J. Murdock, Averaging Methods in Nonlinear Dynamical Systems, Second edition, Applied Mathematical Sciences 59, Springer, New York, 2007. [22] F. Verhulst, Nonlinear Differential Equations and Dynamical Systems, Universitex, Springer, 1991. [23] G.N. Watson, A treatise on the theory of Bessel functions, 2nd ed. Cambridge: Cambridge University Press, 2006.