This is a preprint of: “Limit cycles bifurcating of Kolmogorov systems in R2 and in R3 ”, Jaume Llibre, Yohanna Paulina Martinez Mancilla, Clàudia Valls, Commun. Nonlinear Sci. Numer. Simul., vol. 91, 105401:1–10, 2020. DOI: [10.1016/j.cnsns.2020.10540] LIMIT CYCLES BIFURCATING OF KOLMOGOROV SYSTEMS IN R2 AND IN R3 JAUME LLIBRE1 , Y. PAULINA MARTÍNEZ2,3 AND CLAUDIA VALLS4 Abstract. In this work we consider the Kolmogorov system of degree 3 in R2 and R3 having an equilibrium point in the positive quadrant and octant, respectively. We provide sufficient conditions in order that the equilibrium point will be a Hopf point for the planar case and a zero-Hopf point for the spatial one. We study the limit cycles bifurcating from these equilibria using averaging theory of second and first order, respectively. We note that the equilibrium point is located in the quadrant or octant where the Kolmogorov systems have biological meaning. 1. Introduction and statement of the main results The Lotka–Volterra systems, which are a class of polynomial differential systems of degree 2 in the plane, were developed independently by Alfred J. Lotka in 1925 [6] and Vito Volterra in 1926 [13], they initially proposed these models for studying the interactions between two species. Kolmogorov [4] in 1936 extended these systems to arbitrary dimension and degree, which are now called Kolmogorov systems. The Lotka–Volterra and Kolmogorov systems have been applied to model differ- ent natural phenomena such as the time evolution of conflicting species in biology (see for more details May [9]), the evolution of competition between three species (studied by May and Leonard [8]), the evolution of electrons, ions and neutral species in plasma physics [7], chemical reactions [3], hydrodynamics [2], economics [12], etc. We want to consider the Kolmogorov systems of degree 3 in the plane (resp. in the space) that has an equilibrium (a, b) in the plane (resp. (a, b, c) in the space) in the interior of the first quadrant (resp. octant). We note that the region of ecological interest in Lotka-Volterra and Kolmogorov systems is indeed the first quadrant (resp. octant). It is easy to see that we can consider (a, b) = (1, 1) (resp. (a, b, c) = (1, 1, 1) in the space) doing the rescaling (x, y) → (x/a, y/b) (resp. (x, y, z) → (x/a, y/b, z/c)). In short, in this work we consider the following Kolmogorov systems of degree three in the plane, that is, (1) ẋ = −x(a1 (x − 1) + a2 (y − 1) + a4 (x − 1)2 + a5 (x − 1)(y − 1) + a7 (y − 1)2 ), ẏ = −y(b1 (x − 1) + b2 (y − 1) + b4 (x − 1)2 + b5 (x − 1)(y − 1) + b7 (y − 1)2 ), 2010 Mathematics Subject Classification. Primary 37C10, Secondary 34C05. Key words and phrases. Lotka–Volterra system, Kolmogorov systems, phase portraits, Poincaré ball, Hopf bifurcation, zero-Hopf bifurcation, limit cycle. 1 2 J. LLIBRE, Y.P. MARTÍNEZ AND C. VALLS and in the space, i.e. (2) ẋ = −x(a1 (x − 1) + a2 (y − 1) + a3 (z − 1) + a4 (x − 1)2 + a5 (x − 1)(y − 1)+ a6 (x − 1)(z − 1) + a7 (y − 1)2 + a8 (y − 1)(z − 1) + a9 (z − 1)2 ), ẏ = −y(b1 (x − 1) + b2 (y − 1) + b3 (z − 1) + b4 (x − 1)2 + b5 (x − 1)(y − 1)+ b6 (x − 1)(z − 1) + b7 (y − 1)2 + b8 (y − 1)(z − 1) + b9 (z − 1)2 ), ż = −z(c1 (x − 1) + c2 (y − 1) + c3 (z − 1) + c4 (x − 1)2 + c5 (x − 1)(y − 1)+ c6 (x − 1)(z − 1) + c7 (y − 1)2 + c8 (y − 1)(z − 1) + c9 (z − 1)2 ), where the dot denotes derivative with respect to the time t. As far as we know there are no rigorous analytic studies on the existence of periodic solutions for the Kolmogorov systems (1) and (2) coming from a Hopf bi- furcation for systems (1), and from a zero-Hopf bifurcation for systems (2). In this paper using the averaging theory of first and second order we shall prove the exis- tence of limit cycles for these systems bifurcating from eihter a Hopf or a zero-Hopf equilibrium point. We recall that a Hopf equilibrium point of an autonomous system in R2 is an isolated equilibrium point with a pair of purely imaginary eigenvalues ±iω, ω ∈ R+ . On the other hand, a zero-Hopf equilibrium point of an autonomous system in R3 is an isolated equilibrium point with linear part having one zero eigen- value and a pair of purely imaginary eigenvalues ±iω, ω ∈ R+ . These points are important because under some assumptions a small-amplitude limit cycle emerges from them. We will focus on limit cycles emerging from an equilibrium point (ei- ther Hopf or zero-Hopf) in the interior of the positive quadrant or octant of the Kolmogorov systems (1) and (2), respectively. A limit cycle of systems (1) or (2) is a periodic orbit in the set of all periodic orbits of these differential systems. We say that a limit cycle is unstable if its associated Poincaré map has at least one positive characterize exponent. On the other hand, we say that a limit cycle is stable if its associated Poincaré map has all the characteristic exponents being negative. Our main objective is to provide explicit sufficient conditions for the existence of limit cycles bifurcating from a Hopf equilibrium of a Kolmogorov system of degree 3 in dimension two, and from a zero-Hopf equilibrium of a Kolmogorov system of degree 3 in dimension three. Limit cycles are important in the dynamics of a differential system see for instance [10], and the Hopf bifurcations are also relevant for understanding such dynamics see for example [11]. Let b22 B̂ (3) A= and B = , 2ω 8a32 ω 3 where b22 will be introduced in equation (6), and B̂ = (a21 + a2 (2a3 + b4 ) + a1 (a2 − a4 − 2b5 ))(a31 a5 + a32 b3 + a1 a22 (a3 − b4 ) + a21 a2 (−a4 + b5 )) + (2a31 a5 + a22 (2a2 a3 − a3 a4 + b4 b5 ) + a1 a2 (a24 + 2a3 a5 + a5 b4 + a4 b5 − 2b25 + a2 (a3 − a4 − b4 + b5 )) + a21 (−2a5 (a4 + 2b5 ) + a2 (−a4 + a5 + 3b5 )))ω 2 + (a1 a5 + 2a2 b5 − a5 (a4 + 2b5 ))ω 4 . Our first main result concerning the Kolmogorov systems (1) is the following. Theorem 1. Using the averaging theory of second order, if b1 = −(a21 + ω 2 )/a2 with a2 6= 0, b2 = −a1 and AB < 0, then the Kolmogorov system (1) has one small LIMIT CYCLES BIFURCATING OF KOLMOGOROV SYSTEMS 3 limit cycle bifurcating from the Hopf equilibrium point (1, 1). The limit cycle is stable if A > 0 and unstable if A < 0. Theorem 1 is proved in section 2. The following system in R2 satisfies the hypotheses of Theorem 1. Example 2. The Kolmogorov system ẋ = x(−1 + (x − 1)2 + (x − 1)(y − 1) + (y − 1)2 + y), (4) ẏ = y(−(x − 1) + 2(x − 1)2 + (x − 1)(y − 1) − 10(y − 1)2 ), has one small unstable limit cycle bifurcating from the equilibrium point (1, 1). The details of the example are given in section 2. Our second main result concerns with Kolmogorov systems (2). To state it we need some extra notation. Let α = a21 a3 b1 −a31 b3 +a1 (a3 b1 b2 −b3 (2a2 b1 +ω 2 ))+b1 (a2 (a3 b1 −b2 b3 )+a3 (b22 +ω 2 )), β = −a21 a2 b3 − a22 b1 b3 − a1 a2 (−a3 b1 + b2 b3 ) + a3 b2 (b22 + ω 2 ) − a2 (−2a3 b1 b2 + b3 (b22 + ω 2 )), γ = −a23 b1 + a3 (a1 − b2 )b3 + a2 b23 , δ = −ω(a23 b22 (b22 +ω 2 )+2a2 a3 b2 (a3 b1 b2 −b3 (a1 b2 +b22 +ω 2 ))+a22 (a23 b21 −2a3 b1 (a1 + b2 )b3 + b23 (a21 + 2a1 b2 + b22 + ω 2 ))), S1 = (CS3 )/E, S2 = 2C + c31 /ω, S3 = DS2 + CF . We will also introduce the notation C, D, E, F that due to their complexity will be provided in Appendix B. Theorem 3. Using averaging theory of first order, if c1 = α/γ, c2 = β/γ, c3 = −a1 − b2 , F S2 6= 0,S1 < 0 and CD 6= 0, then the Kolmogorov system (2) has two small limit cycles bifurcating from the zero-Hopf equilibrium point (1, 1, 1). More- over the following holds. (a) If S2 > 0, S3 /F > 0 and DE < 0, or S2 < 0, S3 /F > 0, DE > 0 and (CF + S3 )/D < 0, or S2 > 0, S3 /F < 0, DE > 0 and (CF + S3 )/D < 0, the two limit cycles are unstable. One of them has an unstable manifold of dimension two and the other has a stable manifold and an unstable manifold both of dimension one. (b) If S2 < 0, S3 /F > 0 and DE < 0, or S2 > 0, S3 /F < 0 and DE < 0, the two limit cycles are unstable and both have a stable manifold and an unstable manifold of dimension one. (c) If S2 < 0, S3 /F < 0 and DE < 0, or S2 < 0, S3 /F > 0, DE > 0 and (CF + S3 )/D > 0, or S2 > 0, S3 /F < 0, DE > 0 and (CF + S3 )/D > 0, one limit cycle is stable with a stable manifold of dimension two and the other limit cycle is unstable with a stable manifold and an unstable manifold of dimension one. 4 J. LLIBRE, Y.P. MARTÍNEZ AND C. VALLS (d) If S2 > 0, S3 /F > 0, DE > 0 and (CF + S3 )/D < 0 both limit cycles are unstable with unstable manifolds of dimension two. (e) If S2 < 0, S3 /F < 0, DE > 0 and (CF + S3 )/D > 0 both limit cycles are stable with stable manifolds of dimension two. Theorem 3 is proved in section 3. The following system in R3 satisfies the hypotheses of statement (a) of Theorem 3. Example 4. The Kolmogorov system (5) ẋ = x((y − 1) + (z − 1) + 207/10(x − 1)2 + 17/10(x − 1)(y − 1) − 53/5(y − 1)2 − 29/10(x − 1)(z − 1) − 189/10(y − 1)(z − 1) − 43/10(z − 1)2 ), ẏ = y((x − 1) − 2/5(x − 1)2 + 42/5(x − 1)(y − 1) − 1/2(y − 1)2 + 82/5(x − 1)(z − 1) + 16/5(y − 1)(z − 1) + 3(z − 1)2 ), ż = z(c1 (x − 1) + c2 (y − 1) + c3 (z − 1) − (3/2)(x − 1)2 + 13(x − 1)(y − 1)− 27/10(y − 1)2 + 17/10(x − 1)(z − 1) − 69/5(y − 1)(z − 1) − 87/10(z − 1)2 ), has two small limit cycles bifurcating from the equilibrium point (1, 1, 1). One limit cycle is unstable with an unstable manifold of dimension two and the other is un- stable with a stable manifold and an unstable manifold both of dimension one. The details of the example are given in section 3. See the appendix for a summary of the results on the averaging theory that we shall need for proving Theorems 1 and 3. 2. Proof of Theorem 1 Considering system (1) with a finite equilibria (1, 1) in the interior of the first quadrant. First we impose conditions to the parameters so that the Kolmogorov system (1) has a Hopf point at (1, 1), that is, the Jacobian matrix at (1, 1) has complex con- jugate eigenvalues with zero real part ±iω, ω ∈ R+ . Doing so we get a21 + ω 2 b1 = − , b2 = −a1 , a2 with ω > 0, a2 = 0, and the linearization around the point (1, 1) is   a1 a2 M= −(a21 + ω 2 )/a2 −a1 . We proceed to study the number of limit cycles bifurcating from the Hopf equilib- rium point. To do so, we write b1 and b2 in the form a21 + ω 2 (6) b1 = − + εb11 + ε2 b12 , b2 = −a1 + εb21 + ε2 b22 a2 where ε is a parameter to be taken sufficiently small. We translate the equilibrium point (1, 1) to the origin through the change of variables (x, y) → (x + 1, y + 1) = (X, Y ) and we write the Kolmogorov system into LIMIT CYCLES BIFURCATING OF KOLMOGOROV SYSTEMS 5 its Jordan normal form J using the auxiliary matrix   1 0 P = , such that J = P M P −1 . −a1 /ω −a2 /ω and the variables (u, v) = P (X, Y ), i.e. (7)   a4 u(a1 u + ωv) a5 u̇ = (1 + u) a3 u2 − ωv − + 2 (a1 u + ωv)2 , a2 a2  −1 (a4 u(a1 u + ωv) a5 v̇ = a1 (1 + u)(a3 u2 − ωv − + 2 (a1 u + ωv)2 )+ ω a2 a2 +b5 (a1 u + ωv)2 − a2 (ω(b4 uv + b21 εv + b22 ε2 v + ωu)+  1 2 a1 (b4 u2 + b21 εu + b22 ε2 u − ωv)) + (a2 − a1 u − ωv)(a2 u(b 3 u + ε(b 11 + b 12 ε)) . a22 Now we write system (7) in such a way that conditions of Theorem 5 are satisfied. First we write system (7) in polar coordinates (u, v) = (r cos θ, r sin θ) and get the system of equations (8)  2 ṙ = cos(θ)(r cos(θ) + 1) a5 (a1 r cos(θ) + rω  sin(θ)) − a2 a4 r2 cos2 (θ)(a1 + 2 2 2 2 2 ω tan(θ)) + a2 a3 r cos (θ) − a2 rω sin(θ) /a2 − sin(θ) [a1 (r cos(θ) + 1)(−a2 a4 r2 cos(θ)(a1 cos(θ) + ω sin(θ)) + a5 (a1 r cos(θ) + rω sin(θ))2 + a22 r(a3 r cos2 (θ) − ω sin(θ))) − (a2 − r(a1 cos(θ) + ω sin(θ)))(a1 a2 r(b4 r cos2 (θ)− ω sin(θ)) − b5 (a1 r cos(θ) + rω sin(θ))2 − a22 b3 r2 cos2 (θ) + a2 rω cos(θ)(b4 r sin(θ)+ ω)) + a2 r(a2 − r(a1 cos(θ) + ω sin(θ)))(cos(θ)(a2 b11 − a1 b21 ) − b21 ω sin(θ))+ a2 r2 (a2 − r(a1 cos(θ) + ω sin(θ)))(cos(θ)(a2 b12 − a1 b22 ) − b22 ω sin(θ))]/(a22 ω 2 ) θ̇ = cos2 (θ)[Rω sin(θ)(a21 (a2 − 3a5 ) + a1 a2 (a2 + 2a4 − 2b21  − 2b22 2 − 2b5 )+ a2 (a2 (−a3 + (b11 + b12 ) + b4 ) − ω 2 )) + a22 (a1 (b21 + b22 ) − a2 (b11 + b12 )+ ω 2 ) + R2 ω 2 2 sin2 (θ)(3a1 (b5 − a5 ) + a2 (a4 − b4 ))) − a1 R2 2 cos4 (θ)(a21 (a5 − b5 )+ a1 a2 (b4 − a4 ) + a22 (a3 − b3 )) − R cos3 (θ)(a31 a5 + Rω sin(θ)(3a21 (a5 − b5 )+ 2a1 a2 (b4 − a4 ) + a22 (a3 − b3 )) + a21 a2 (−a4 + (b21 + b22 ) + b5 ) + a1 a2 (a2 (a3 − (b11 + b12 ) − b4 ) + ω 2 ) + a32 b3 ) + 41 ω sin(2θ)(2Rω sin(θ)(a1 (a2 − 3a5 ) + a2 (a2 + a4 − b21  − b22 2 − b5 )) + 2a22 b21 + 2a22 b22  + R2 ω 2 (a5 − b5 ) cos(2θ)− a5 R2 ω 2  + b5 R2 ω 2 ) + ω 2 sin2 (θ)(a22 − a5 Rω sin(θ)]/(a22 ω) Note that system (8) is periodic in the variable θ with period 2π. We write system (8) by taking as the new independent variable the angular variable θ and we obtain r0 = θ̇ṙ where the prime denotes derivative with respect to θ. We now consider the new variable R = r/ε. We compute R0 and we develop the new equation for R0 in power series in the variable ε up to second order in the form (9) R0 = εF1 + ε2 F2 + O(ε3 ), where F1 = −R[−(a22 a3 −a1 a2 a4 +a21 a5 )Rω cos3 θ +R(a31 a5 +a21 a2 (−a4 +b5 )+a2 (a22 b3 + a2 ω 2 +a4 ω 2 )+a1 (a22 (a3 −b4 )+a2 ω 2 −2a5 ω 2 )) cos2 θ sin θ+ω sin2 θ(−a22 b21 −(a1 a2 − a1 a5 − a2 b5 )Rω sin θ) − cos θ sin θ(a22 (−a2 b11 + a1 b21 ) + Rω(a21 (a2 − 2a5 ) + a22 b4 + a1 a2 (a2 + a4 − 2b5 ) − a2 ω 2 + a5 ω 2 ) sin θ)]/(a22 ω 2 ), 6 J. LLIBRE, Y.P. MARTÍNEZ AND C. VALLS and F2 = R[(a22 a3 − a1 a2 a4 + a21 a5 )Rω cos3 θ − R(a31 a5 + a21 a2 (−a4 + b5 ) + a2 (a22 b3 + a2 ω 2 + a4 ω 2 ) + a1 (a22 (a3 − b4 ) + a2 ω 2 − 2a5 ω 2 )) cos2 θ sin θ + ω sin θ2 (a22 b21 + (a1 a2 − a1 a5 − a2 b5 )Rω sin θ) + cos θ sin θ(a22 (−a2 b11 + a1 b21 ) + Rω(a21 (a2 − 2a5 ) + a22 b4 + a1 a2 (a2 +a4 −2b5 )−a2 ω 2 +a5 ω 2 ) sin θ))(R(a31 a5 +a32 b3 +a21 a2 (−a4 +b5 )+a1 a2 (a2 a3 − a2 b4 + ω 2 )) cos3 θ + a5 Rω 3 sin3 θ − ω cos θ sin θ(a22 b21 + (a1 (a2 − 3a5 ) + a2 (a2 + a4 − b5 ))Rω sin θ) + cos2 θ(a22 (a2 b11 − a1 b21 ) + Rω(−a21 (a2 − 3a5 ) − a1 a2 (a2 + 2a4 − 2b5 ) + a2 (a2 a3 − a2 b4 + ω 2 )) sin θ)) + a22 ((a22 a3 − a1 a2 a4 + a21 a5 )R2 ω cos4 θ − R2 (a21 a2 (−a4 + b4 )+a31 (a5 −b5 )+a2 a4 ω 2 +a1 (a22 (a3 −b3 )−2a5 ω 2 )) cos3 θ sin θ−R cos2 θ sin θ(a1 a2 (−a2 b1 1+ a1 b2 1) − Rω(a22 b3 + a1 a2 (a4 − 2b4 ) + a21 (−2a5 + 3b5 ) + a5 ω 2 ) sin θ) + ω sin2 θ(a22 b22 − a2 b21 Rω sin θ+b5 R2 ω 2 sin2 θ)−cos θ sin θ(a22 (a2 b12 −a1 b22 )−a2 (a2 b11 −2a1 b21 )Rω sin θ+ (a1 a5 + a2 b4 − 3a1 b5 )R2 ω 2 sin2 θ))]/(a42 ω 4 ). Note that the differential system (9) is in normal form (12) for applying the averaging theory with T = 2π, x = R, t = θ and ε3 R(θ, x, ε) = O(ε3 ). We also observe that F1 is C 2 in x and 2π-periodic in θ in an open interval I = {R : 0 < R < R̄} for some R̄ > 0. Applying Theorem 5 (see Appendix A) we obtain the averaging function of first order b21 f1 (R) = R. 2ω We can see that f1 (R) has no solution in I. Therefore the averaging method of first order does not provide any small limit cycle bifurcating from the origin. We set b21 = 0 such that f1 (R) ≡ 0 and then can apply the averaging theory of second order. After some calculations using F2 and that  y1 = − 12aR2 ω2 −Rω(9 sin(θ) + sin(3θ)) a21 a5 − a1 a2 a4 + a22 a3 − 2 2 sin2 (θ)(2Rω sin(θ)(a21 (a2 − 2a5 ) + a1 a2 (a2 + a4 − 2b5 ) + a22 b4 − a2 ω 2 +a5 ω 2 ) − 3a32 b11 ) − 4R(cos3 (θ) − 1)(a31 a5 + a21 a2 (b5 − a4 ) + a1 a22 (a3 − b4 )+ a1 ω 2 (a2 − 2a5 ) + a32 b3 + a2 ω 2 (a2 + a4 )) + 16Rω 2 sin4 ( θ2 )(cos(θ) + 2)(a2 b5 + a1 (a5 − a2 ))) , following (13) we get that the averaging function of second order is f2 (R) = R(A + B R2 ), where p A and B are given in (3). Thus, f2 (R) has one positive real root R∗ = ∗ + −A/B in I if 0 < −A/B and in this case it holds that df2 /dR(R ) = −2A 6= 0. It follows then from Theorem 5 that for |ε| sufficiently p small, system (9) has one 2π-periodic solution R∗ (θ, ε) such that R∗ (θ, ε) → −A/B when ε → 0. Moreover, this periodic solution is stable if A > 0 and unstable if A < 0. Now we shall go back through the changes of variables and we obtain a periodic solution (x(t, ε), y(t, ε)) bifurcating from (1, 1) with a period tending to 2π when ε → 0. Moreover, (x(t, ε), y(t, ε)) = (1, 1) + O(ε). This completes the proof of the theorem. 2.1. Details of Example 2. Consider system (4) with a1 = 1, a2 = 0, b1 = −1 and b2 = 0. Note that it is in the form (1) and the eigenvalues associated to the Jacobian matrix around (1, 1) are equal to ±i. So, (1, 1) is a Hopf equilibrium point. LIMIT CYCLES BIFURCATING OF KOLMOGOROV SYSTEMS 7 For this we perturb the parameters b1 and b2 as in (6) taking b1 = −1 + b11 ε + b12 ε2 and b2 = b22 ε2 , for some parameter ε to be taken sufficiently small. After translating the Hopf equilibrium point to the origin, introducing polar coordinates and setting r = εR we get  Ṙ = 12 R sin θ(Rε sin θ + 1) 2 cos θ b11 ε + b12 ε2 − 1 + 2b22 ε2 sin θ + Rε sin(2θ)+ 12Rε cos(2θ) − 8Rε) + 2 cos θ(Rε cos θ + 1)(sin θ + Rε sin θ cos θ + Rε)) θ̇ = 12 cos θ(Rε sin θ + 1) 2 cos θ b11 ε + b12 ε2 − 1 + 2b22 ε2 sin θ + Rε sin(2θ)+ 12Rε cos(2θ) − 8Rε) − 2 sin θ(Rε cos θ + 1)(sin θ + Rε sin θ cos θ + Rε)) . Finally, we write R0 = Ṙ/θ̇ as in (9) as R0 = −Rε(cos θ(b11 sin θ + R) + R sin θ(7 cos(2θ) − 3)) − 81 Rε2 2b211 sin(2θ)+ +b211 sin(4θ) − 16b11 R sin θ − 2b11 R sin(3θ) + 14b11 R sin(5θ)+  +6b11 R cos θ + 2b11 R cos(3θ) + 4b12 sin(2θ) − 4 b22 − 3R2 cos(2θ) +  4b22 + 109R2 sin(2θ) − 112R2 sin(4θ) + 49R2 sin(6θ) + 16R2 cos(4θ) − 4R2 . The averaging function of first order is identically zero. The averaging function of second order is f2 (R) = −R(b22 − R2 )/2. √ Then for b22 > 0 we have the solution R∗ = b22 6= 0 and so df2 /dR = b22 > 0. Therefore, there exists a unique unstable limit cycle bifurcating from the Hopf-point (1, 1). 3. Proof of Theorem 3 The proof of Theorem 3 is analogous to the proof of Theorem 1 in section 2, so we will follow the structure of the previous section. Consider system (2) with (1, 1, 1) a finite equilibrium point. The conditions to have a zero-Hopf equilibrium point at (1, 1, 1) are c1 = α/γ, c2 = β/γ, c3 = −a1 − b2 . The Jacobian matrix evaluated at (1, 1, 1) has the eigenvalues 0, ±iω, ω ∈ R+ and so it is a zero-Hopf equilibrium point. To study the number of limit cycles bifurcating from this equilibrium point we take c1 = α/γ + εc11 , c2 = β/γ + εc21 , c3 = −a1 − b2 + εc31 , being ε a small parameter to be taken sufficiently small. Doing the same steps detailed in the proof of Theorem 1 for the Kolmogorov system (1) in R2 we arrive to the system in the variables (u, v, w) such that the zero-Hopf equilibrium point is located at the origin and its linear part is in Jordan normal form. Now we want to write the system in such a way that conditions of Theorem 5 are satisfied. 8 J. LLIBRE, Y.P. MARTÍNEZ AND C. VALLS For this we introduce cylindrical coordinates (u, v, w) → (r sin θ, r cos θ, w) and get the associated vector field U = (ṙ(r, θ, w), θ̇(r, θ, w), ẇ(r, θ, w)). Now we intro- duce the new variables (r, w) → (Rε, W ε) and we obtain a system in the form Ṙ = ṙ/ε = G11 (R, θ, W )ε + O(ε2 ), θ̇ = ω + O(ε), Ẇ = ẇ/ε = G12 (R, θ, W ) + O(ε2 ) where Gij are 2π-periodic in the variable θ of class C 2 and O(ε2 ) is also of class C 2 and 2π-periodic in the variable θ. We take θ as the new independent variable and write R0 = Ṙ/θ̇, W 0 = Ẇ /θ̇. Doing so, we get (10) R0 = εF11 + O(ε2 ), W 0 = εF12 + O(ε2 ). Note that the differential system (10) is in normal form (12) for applying the averaging theory with T = 2π, x = (R, W ), t = θ and εR(θ, x, ε) = O(ε2 ). We also remark that F11 , F12 and R are C 2 in x and 2π-periodic in θ in an open interval I × R. Applying Theorem 5 we obtain the averaging function of first order has the form f = (f11 (R, W ), f12 (R, W )) where Z 2π 1 R̂ f11 (R, W ) = dθ = R(C + DW ) 2π 0 ω and Z 2π 1 Ŵ f12 (R, W ) = dθ = ER2 − S2 W + F W 2 , 2π 0 ω with C, D, E and F explicit in Appendix B. We can see that f11 (R, W ) = 0 has the solutions R1 = 0 and W2 = −C/D if CD 6= 0. For R1 = 0 the function f12 (R, W ) = 0 has the solution W1 = S2 /F if S2 F 6=√ 0 and for W = W2 the function f12 (R, W ) = 0 has the positive solution R2 = −S1 /D if S1 < 0. Now we compute the Jacobian matrix of f and we get   C + DW DR . 2ER −S2 + 2F W By evaluating the determinant of the Jacobian matrix at the solution (R1 , W1 ) we get that it is equal to S2 S3 /F and at the solution (R2 , W2 ) we get that it is equal to −2CS3 /D. It follows from Theorem 5 that if S3 6= 0 then for |ε| sufficiently small, system (10) has two 2π-periodic solutions (R1 (θ, ε), W1 (θ, ε)) and (R2 (θ, ε), W2 (θ, ε)) such that (Rj (θ, ε), Wj (θ, ε)) → (Rj , Wj ) for j = 1, 2 when ε → 0. Moreover, the Jacobian matrix evaluated at the periodic solution (R1 , W1 ) has eigenvalues equal to S2 and S3 /F . Since the eigenvalues of the Jacobian matrix evaluated at the periodic solutions (R1 , W1 ) provide the stability of the fixed point corresponding to the Poincaré map defined in a neighborhood of the periodic solu- tion, if S2 > 0 and S3 /F > 0 then the periodic solution is unstable with an unstable manifold of dimension two. If S2 < 0 and S3 /F > 0, or S2 > 0 and S3 /F < 0 then the periodic solution is unstable with a stable and an unstable manifold of dimension one. Finally, if S2 < 0 and S3 /F < 0 then the periodic solution is stable with a stable manifold of dimension two. LIMIT CYCLES BIFURCATING OF KOLMOGOROV SYSTEMS 9 On the other hand,√the Jacobian matrix evaluated at (R2 , W2 ) has eigenvalues CF + S3 ± S4 equal to − where S4 = (CF +S3 )2 +8CDS3 . Since S1 = CS3 /E < 2D 0, then taking into account that CDS3 = DES1 we get that if DE < 0 then the periodic solution is unstable with a stable and an unstable manifold of dimension one. If DE > 0 and (CF + S3 )/D > 0 then the periodic solution is stable with a stable manifold of dimension two. Finally, if DE > 0 and (CF + S3 )/D < 0 then the periodic solution is unstable with a unstable manifold of dimension two. Combining the above information of the eigenvalues of the Jacobian matrix for both (R1 , W1 ) and (R2 , W2 ) we get statements (a)–(e) in the theorem. Now we shall go back through the changes of variables and we obtain two periodic solutions (xj (t, ε), y(t, ε), zj (t, ε)) (j = 1, 2) bifurcating from (1, 1, 1) with a period tending to 2π when ε → 0. Moreover, (xj (t, ε), y(t, ε), zj (t, ε)) = (1, 1, 1) + O(ε) for j = 1, 2. This completes the proof of the theorem. 3.1. Details on Example 4. . Take system (5). It is in the form (2), with c1 = −2, c2 = 0 and c3 = 0. Moreover, the singular point (1, 1, 1) is a zero-Hopf equilibrium point (the eigenvalues of the Jacobian matrix at this point are 0, ±i). We write c1 = −2 + c11 ε, c2 = c21 ε, c3 = c31 ε = −ε/5, with ε a small parameter to be taken sufficiently small. Thus, following the pro- cedure we did it in the proof of Theorem 3, translating the zero-Hopf equilibrium point at the origin, writing the Jacobian matrix at this point in Jordan normal form, introducing cylinder coordinates and taking the variables R, W applying the averaging theory of first order we get the functions (11) f11 (R, W ) = R(1 + 2W ), f12 (R, W ) = (−R2 + 2/5W (−9 + 10W ))/2. p From (11) we have two solutions namely s1 = (0, 9/10) and s2 = ( 14/5, −1/2), whose determinant of the Jacobian matrix is equal to 126/25 and 28/5, respectively. Moreover, the eigenvalues of s1 are 14/5 and 9/5 so it is unstable with an unstable manifold √ of dimension two. On √ the other hand, the eigenvalues of s2 are −(38 + 2 921)/20 > 0 and −(38 − 2 921)/20 < 0 so s2 is unstable with a stable and an unstable manifolds of dimension one. Then, through the averaging method of first order we have detected two small limit cycles bifurcating from the zero-Hopf equilibrium point (1, 1, 1) of system (5) when ε = 0. Appendix A. Averaging Theory We summarize the averaging theory of second order which provides sufficient conditions for the existence of periodic orbits for a periodic differential system depending on small parameters. See [1] for additional details and of the prof of the result stated in this appendix. Theorem 5. Consider the differential system (12) ẋ(t) = εF1 (t, x) + ε2 F2 (t, x) + ε3 R(t, x, ε), 10 J. LLIBRE, Y.P. MARTÍNEZ AND C. VALLS where F1 , F2 : R × D → Rn , R : R × D × (−εf , εf ) → Rn are continuous and T - periodic functions in the first variable, and D is an open subset of Rn . Assume that the following assumptions hold: (a) F1 (t, ·) ∈ C 1 (D) for all t ∈ R, F1 , F2 , R and Dx F1 are locally Lipschitz with respect to x, and R is differentiable with respect to ε. We define f1 , f2 : D → Rn as Z T f1 (z) = F1 (s, z) ds, 0 (13) Z T f2 (z) = [Dz F1 (s, z)y1 (s, z) + F2 (s, z)] ds 0 where Z s y1 (s, z) = F1 (t, z) dt. 0 (b) For V ⊂ D an open and bounded set and for each ε ∈ (−εf , εf ) \ {0}, there exists a ∈ V such that if f1 (z) 6≡ 0 then f1 (a) = 0 and dB (f, a) 6= 0, where dB (f, a) is the Brouwer degree of the function f1 : V → Rn at the fixed point a, and if f1 (z) ≡ 0 and f2 (z) 6≡ 0, then f2 (a) = 0 and dB (f, a) 6= 0. Then for |ε| > 0 sufficiently small, there exists a T -periodic solution φ(t, ε) of the system such that φ(0, ε) → a when ε → 0. The kind of stability or instability of the limit cycle φ(t, ε) is given by the eigenvalues of the Jacobian matrix D2 (f1 (z) + εf2 (z))|z=a . Note that a sufficient condition for the showing that the Brouwer degree of a function f at a fixed point a is nonzero, is that the Jacobian of the function f at a (when it is defined, it is not zero), see [5]. When the first part of assumption (b) holds Theorem 5 provides the averaging theory of first order and it provides the averaging theory of second order when the second part of assumption (b) holds. Appendix B. notations In this section we introduce the notation introduced in the statement and proof of Theorem 3. C = (a3 b2 c11 − a2 b3 c11 − a3 b1 c21 + a1 b3 c21 + a2 b1 c31 − a1 b2 c31 + c31 ω 2 )/(2ω 3 ), D = (2a41 b3 (a9 b22 +b3 (−a8 b2 +a7 b3 ))−2a43 b1 (b22 c4 −b1 b2 c5 +b21 c7 )+a22 b3 (2a2 (a9 b31 − a6 b21 b3 + a4 b1 b23 + b2 b23 b4 − b1 b2 b3 b6 + b21 b2 b9 + b33 c4 − b1 b23 c6 + b21 b3 c9 ) + (2a9 b21 − b3 (a6 b1 + b3 (b5 + c6 ) + b1 (b2 − b8 − 2c9 )))ω 2 ) + a33 (−2a7 b31 b2 − 2a4 b1 b32 − 2a2 b1 b22 b4 − 2b42 b4 + 2a2 b21 b2 b5 + 2b1 b32 b5 − 2a2 b31 b7 − 2b21 b22 b7 + 4a2 b1 b2 b3 c4 − 2b32 b3 c4 − 2a2 b21 b3 c5 + 2b1 b22 b3 c5 −2a2 b21 b2 c6 −2b21 b2 b3 c7 +2a2 b31 c8 −2a4 b1 b2 ω 2 +b21 b2 ω 2 −2b22 b4 ω 2 +b1 b2 b5 ω 2 − b1 b2 c6 ω 2 + b21 c8 ω 2 + a5 b21 (2b22 + ω 2 )) − a2 a3 (2a22 b1 (b23 b4 − b1 b3 b6 + b21 b9 ) + b3 (2a8 b21 − a6 b1 b2 − b1 b22 − a5 b1 b3 + b1 b2 b3 − 2b2 b3 b5 + 2b1 b3 b7 + b1 b2 b8 − 2b2 b3 c6 + b1 b3 c8 + 2b1 b2 c9 )ω 2 + 2a2 (a9 b31 b2 + a8 b31 b3 − 2a6 b21 b2 b3 − a5 b21 b23 + 3a4 b1 b2 b23 + 3b22 b23 b4 − b1 b2 b23 b5 − 2b1 b22 b3 b6 + b21 b2 b3 b8 + b21 b22 b9 + 3b2 b33 c4 − b1 b33 c5 − 2b1 b2 b23 c6 + b21 b23 c8 + b21 b2 b3 c9 +b23 b4 ω 2 −b1 b3 b6 ω 2 +b21 b9 ω 2 ))+a23 (2a22 b1 (b2 (2b3 b4 −b1 b6 )−b23 c4 +b1 b3 (−b5 + c6 ) + b21 (b8 − c9 )) + b3 (2a7 b21 + b2 (−a5 b1 − b2 (b5 + c6 ) + b1 (b2 + 2b7 + c8 )))ω 2 + LIMIT CYCLES BIFURCATING OF KOLMOGOROV SYSTEMS 11 a2 (2a8 b31 b2 + 2a7 b31 b3 − 4a5 b21 b2 b3 + 6a4 b1 b22 b3 + 6b32 b3 b4 − 4b1 b22 b3 b5 − 2b1 b32 b6 + 2b21 b2 b3 b7 +2b21 b22 b8 +6b22 b23 c4 −4b1 b2 b23 c5 −2b1 b22 b3 c6 +2b21 b23 c7 +2b21 b2 b3 c8 +b21 b2 ω 2 + 2a4 b1 b3 ω 2 + 4b2 b3 b4 ω 2 − b1 b3 b5 ω 2 − 2b1 b2 b6 ω 2 + b21 b8 ω 2 + b1 b3 c6 ω 2 − 2b21 c9 ω 2 − a6 b21 (2b22 + ω 2 ))) + a31 (2a2 b3 (−2a9 b1 b2 + a8 b1 b3 + a6 b2 b3 − a5 b23 + b23 b7 − b2 b3 b8 + b22 b9 ) + a3 (−2a9 b1 b22 + b3 (4a8 b1 b2 − 2a6 b22 − 6a7 b1 b3 + 2a5 b2 b3 + 2b23 c7 − 2b2 b3 c8 + 2b22 c9 + b2 ω 2 ))) + a21 (a23 (−2a8 b21 b2 + 2a6 b1 b22 + 6a7 b21 b3 − 4a5 b1 b2 b3 + 2a4 b22 b3 + 2b2 b23 c5 − 2b22 b3 c6 − 6b1 b23 c7 + 4b1 b2 b3 c8 − 2b1 b22 c9 − b1 b2 ω 2 + b2 b3 ω 2 ) + b3 (2a22 (a9 b21 − a6 b1 b3 + a4 b23 − b23 b5 + b2 b3 b6 + b1 b3 b8 − 2b1 b2 b9 ) + 2(a9 b22 − a8 b2 b3 + a7 b23 )ω 2 + a2 (2a9 b1 b22 − 2a8 b1 b2 b3 + 2a7 b1 b23 + 2b2 b23 b7 − 2b22 b3 b8 + 2b32 b9 + 2b33 c7 − 2b2 b23 c8 + 2b22 b3 c9 +b2 b3 ω 2 ))+a3 (−2a9 b1 b32 +2a8 b1 b22 b3 −2a7 b1 b2 b23 −2b22 b23 b7 +2b32 b3 b8 −2b42 b9 − 2b2 b33 c7 + 2b22 b23 c8 − 2b32 b3 c9 − a6 b2 b3 ω 2 + a5 b23 ω 2 + b2 b23 ω 2 + b2 b3 b8 ω 2 − 2b22 b9 ω 2 + b23 c8 ω 2 − 2b2 b3 c9 ω 2 + a2 (4a9 b21 b2 − 4a8 b21 b3 + 4a5 b1 b23 − 4a4 b2 b23 + 2b2 b23 b5 − 2b22 b3 b6 − 6b1 b23 b7 + 4b1 b2 b3 b8 − 2b1 b22 b9 − 2b33 c5 + 2b2 b23 c6 + 2b1 b23 c8 − 4b1 b2 b3 c9 − b1 b3 ω 2 − b23 ω 2 ))) − a1 (a33 (2a7 b31 − 2a5 b21 b2 + 2a4 b1 b22 − 2b22 b3 c4 + 4b1 b2 b3 c5 − 2b1 b22 c6 − 6b21 b3 c7 + 2b21 b2 c8 + b1 b2 ω 2 ) + a23 (2a8 b21 b22 − 2a6 b1 b32 − 4a7 b21 b2 b3 + 2a5 b1 b22 b3 + 2b32 b3 b5 − 2b42 b6 − 4b1 b22 b3 b7 + 2b1 b32 b8 + 2b22 b23 c5 − 2b32 b3 c6 − 4b1 b2 b23 c7 + 2b1 b22 b3 c8 − a6 b1 b2 ω 2 + b1 b22 ω 2 + 2a5 b1 b3 ω 2 − 2a4 b2 b3 ω 2 + 2b1 b2 b3 ω 2 + b22 b3 ω 2 + b2 b3 b5 ω 2 − 2b22 b6 ω 2 + b1 b2 b8 ω 2 − b2 b3 c6 ω 2 + 2b1 b3 c8 ω 2 − 2b1 b2 c9 ω 2 − a2 (2a8 b31 − 2a6 b21 b2 − 2a5 b21 b3 + 4a4 b1 b2 b3 + 2b22 b3 b4 − 4b1 b2 b3 b5 + 2b1 b22 b6 + 6b21 b3 b7 − 2b21 b2 b8 − 4b2 b23 c4 + 4b1 b23 c5 − 4b21 b3 c8 + 4b21 b2 c9 + b21 ω 2 + b1 b3 ω 2 )) + a2 b3 (−2a22 (b23 b4 − b1 b3 b6 + b21 b9 ) + (4a9 b1 b2 − b3 (2a8 b1 + a6 b2 + b22 − a5 b3 + b2 b3 + 2b3 b7 − b2 b8 + b3 c8 − 2b2 c9 ))ω 2 + a2 (4a9 b21 b2 − 2a8 b21 b3 − 2a6 b1 b2 b3 + 2a5 b1 b23 + 2b2 b23 b5 − 2b22 b3 b6 − 2b1 b2 b3 b8 + 4b1 b22 b9 + 2b33 c5 − 2b2 b23 c6 − 2b1 b23 c8 + 4b1 b2 b3 c9 + b1 b3 ω 2 + b23 ω 2 )) + a3 (2a22 (a9 b31 − a6 b21 b3 + a4 b1 b23 + 2b2 b23 b4 − 2b1 b23 b5 + 2b21 b3 b8 − 2b21 b2 b9 − b33 c4 + b1 b23 c6 − b21 b3 c9 ) + b3 (−2a8 b1 b2 + a6 b22 + b32 + 4a7 b1 b3 −a5 b2 b3 +b22 b3 +2b2 b3 b7 −b22 b8 +b2 b3 c8 −2b22 c9 )ω 2 +a2 (−4a9 b21 b22 +4a7 b21 b23 − 4a5 b1 b2 b23 − 4b22 b23 b5 + 4b32 b3 b6 + 4b1 b2 b23 b7 − 4b1 b32 b9 − 4b2 b33 c5 + 4b22 b23 c6 + 4b1 b33 c7 − 4b1 b22 b3 c9 +2a4 b23 ω 2 −2b2 b23 ω 2 −b23 b5 ω 2 +2b2 b3 b6 ω 2 +b1 b3 b8 ω 2 −4b1 b2 b9 ω 2 +b23 c6 ω 2 − 2b1 b3 c9 ω 2 + a6 b1 b3 (4b22 − ω 2 )))))/(2(a23 b1 + a3 (−a1 + b2 )b3 − a2 b23 )ω 5 ), E = ((a2 b1 −a1 b2 )2 (a43 b1 (b22 c4 −b1 b2 c5 +b21 c7 +c4 ω 2 )+a33 (a7 b31 b2 −a5 b21 b22 +a4 b1 b32 + a2 b1 b22 b4 +b42 b4 −a2 b21 b2 b5 −b1 b32 b5 +a2 b31 b7 +b21 b22 b7 −2a2 b1 b2 b3 c4 +b32 b3 c4 +a2 b21 b3 c5 − b1 b22 b3 c5 + a2 b21 b2 c6 + b21 b2 b3 c7 − a2 b31 c8 + a4 b1 b2 ω 2 + a2 b1 b4 ω 2 + 2b22 b4 ω 2 − b1 b2 b5 ω 2 + b21 b7 ω 2 +b2 b3 c4 ω 2 +b1 b3 c5 ω 2 −b21 c8 ω 2 +b4 ω 4 )−a41 b3 (b3 (−a8 b2 +a7 b3 )+a9 (b22 +ω 2 ))− a23 (a22 b1 (b2 (2b3 b4 − b1 b6 ) − b23 c4 + b1 b3 (−b5 + c6 ) + b21 (b8 − c9 )) + ω 2 (a8 b21 b2 + a7 b21 b3 − 2a5 b1 b2 b3 + a4 b22 b3 − b22 b3 b5 + b1 b22 b8 − b2 b23 c5 − b1 b23 c7 + 2b1 b2 b3 c8 − b1 b22 c9 + a4 b3 ω 2 − b3 b5 ω 2 + b1 b8 ω 2 − b1 c9 ω 2 ) + a2 (a8 b31 b2 − a6 b21 b22 + a7 b31 b3 − 2a5 b21 b2 b3 + 3a4 b1 b22 b3 + 3b32 b3 b4 − 2b1 b22 b3 b5 − b1 b32 b6 + b21 b2 b3 b7 + b21 b22 b8 + 3b22 b23 c4 − 2b1 b2 b23 c5 − b1 b22 b3 c6 + b21 b23 c7 +b21 b2 b3 c8 +a4 b1 b3 ω 2 +3b2 b3 b4 ω 2 −2b1 b3 b5 ω 2 −b1 b2 b6 ω 2 +2b21 b8 ω 2 +b23 c4 ω 2 + b1 b3 c6 ω 2 − 2b21 c9 ω 2 )) + a3 (a32 b1 (b23 b4 − b1 b3 b6 + b21 b9 ) + a2 ω 2 (2a9 b21 b2 + 2a8 b21 b3 − 2a6 b1 b2 b3 − 2a5 b1 b23 + 2a4 b2 b23 − b2 b23 b5 − b22 b3 b6 + b1 b23 b7 + 3b1 b22 b9 − b33 c5 − b2 b23 c6 + b1 b23 c8 + 2b1 b2 b3 c9 − b3 b6 ω 2 + 3b1 b9 ω 2 ) + a22 (a9 b31 b2 + a8 b31 b3 − 2a6 b21 b2 b3 − a5 b21 b23 + 3a4 b1 b2 b23 + 3b22 b23 b4 − b1 b2 b23 b5 − 2b1 b22 b3 b6 + b21 b2 b3 b8 + b21 b22 b9 + 3b2 b33 c4 − b1 b33 c5 − 2b1 b2 b23 c6 +b21 b23 c8 +b21 b2 b3 c9 +b23 b4 ω 2 −2b1 b3 b6 ω 2 +3b21 b9 ω 2 )+ω 2 (a7 b1 b2 b23 +b22 b23 b7 − b32 b3 b8 + b42 b9 + b2 b33 c7 − b22 b23 c8 + b32 b3 c9 − a5 b23 ω 2 + b23 b7 ω 2 − b2 b3 b8 ω 2 + 2b22 b9 ω 2 + b2 b3 c9 ω 2 +b9 ω 4 +a8 b1 b3 (−b22 +ω 2 )+a9 b1 b2 (b22 +ω 2 )))−b3 (a32 (a9 b31 −a6 b21 b3 +a4 b1 b23 + b2 b23 b4 −b1 b2 b3 b6 +b21 b2 b9 +b33 c4 −b1 b23 c6 +b21 b3 c9 )+a22 (3a9 b21 −2a6 b1 b3 +a4 b23 −b2 b3 b6 + 2b1 b2 b9 − b23 c6 + 2b1 b3 c9 )ω 2 + ω 4 (b3 (−a8 b2 + a7 b3 ) + a9 (b22 + ω 2 )) + a2 ω 2 (−a8 b1 b2 b3 + a7 b1 b23 + b2 b23 b7 − b22 b3 b8 + b32 b9 + b33 c7 − b2 b23 c8 + b22 b3 c9 − a6 b3 ω 2 + b2 b9 ω 2 + b3 c9 ω 2 + a9 b1 (b22 + 3ω 2 ))) + a31 (a2 b3 (2a9 b1 b2 − a8 b1 b3 − a6 b2 b3 + a5 b23 − b23 b7 + b2 b3 b8 − b22 b9 − 12 J. LLIBRE, Y.P. MARTÍNEZ AND C. VALLS b9 ω 2 )+a3 (a9 b1 (b22 +ω 2 )−b3 (2a8 b1 b2 −3a7 b1 b3 +a5 b2 b3 +b23 c7 −b2 b3 c8 +b22 c9 +c9 ω 2 − a6 (b22 +ω 2 ))))+a21 (a23 (a8 b21 b2 −3a7 b21 b3 +2a5 b1 b2 b3 −a4 b22 b3 −b2 b23 c5 +b22 b3 c6 +3b1 b23 c7 − 2b1 b2 b3 c8 + b1 b22 c9 − a4 b3 ω 2 + b3 c6 ω 2 + b1 c9 ω 2 − a6 b1 (b22 + ω 2 )) + a3 (−a8 b1 b22 b3 + a7 b1 b2 b23 + b22 b23 b7 − b32 b3 b8 + b42 b9 + b2 b33 c7 − b22 b23 c8 + b32 b3 c9 + a8 b1 b3 ω 2 − a5 b23 ω 2 + b23 b7 ω 2 − b2 b3 b8 ω 2 + 2b22 b9 ω 2 + b2 b3 c9 ω 2 + b9 ω 4 + a9 b1 b2 (b22 + ω 2 ) + a2 (−2a9 b21 b2 + 2a8 b21 b3 − 2a5 b1 b23 + 2a4 b2 b23 − b2 b23 b5 + b22 b3 b6 + 3b1 b23 b7 − 2b1 b2 b3 b8 + b1 b22 b9 + b33 c5 − b2 b23 c6 − b1 b23 c8 + 2b1 b2 b3 c9 + b3 b6 ω 2 + b1 b9 ω 2 )) − b3 (a22 (a9 b21 − a6 b1 b3 + a4 b23 − b23 b5 + b2 b3 b6 + b1 b3 b8 − 2b1 b2 b9 ) + 2ω 2 (b3 (−a8 b2 + a7 b3 ) + a9 (b22 + ω 2 )) + a2 (−a8 b1 b2 b3 + a7 b1 b23 + b2 b23 b7 − b22 b3 b8 + b32 b9 + b33 c7 − b2 b23 c8 + b22 b3 c9 − a6 b3 ω 2 + b2 b9 ω 2 + b3 c9 ω 2 + a9 b1 (b22 + 3ω 2 )))) + a1 (a33 (a7 b31 − a5 b21 b2 + a4 b1 b22 − b22 b3 c4 + 2b1 b2 b3 c5 − b1 b22 c6 − 3b21 b3 c7 + b21 b2 c8 + a4 b1 ω 2 − b3 c4 ω 2 − b1 c6 ω 2 ) − a23 (a6 b1 b32 + 2a7 b21 b2 b3 − a5 b1 b22 b3 − b32 b3 b5 +b42 b6 +2b1 b22 b3 b7 −b1 b32 b8 −b22 b23 c5 +b32 b3 c6 +2b1 b2 b23 c7 −b1 b22 b3 c8 +a6 b1 b2 ω 2 − a5 b1 b3 ω 2 −b2 b3 b5 ω 2 +2b22 b6 ω 2 +2b1 b3 b7 ω 2 −b1 b2 b8 ω 2 +b23 c5 ω 2 +b2 b3 c6 ω 2 −b1 b3 c8 ω 2 + b6 ω 4 + a8 b21 (−b22 + ω 2 ) + a2 (a8 b31 − a6 b21 b2 − a5 b21 b3 + 2a4 b1 b2 b3 + b22 b3 b4 − 2b1 b2 b3 b5 + b1 b22 b6 +3b21 b3 b7 −b21 b2 b8 −2b2 b23 c4 +2b1 b23 c5 −2b21 b3 c8 +2b21 b2 c9 +b3 b4 ω 2 +b1 b6 ω 2 ))− a2 b3 (a22 (b23 b4 −b1 b3 b6 +b21 b9 )+ω 2 (−2a9 b1 b2 +a8 b1 b3 +a6 b2 b3 −a5 b23 +b23 b7 −b2 b3 b8 + b22 b9 +b9 ω 2 )+a2 (−2a9 b21 b2 +a8 b21 b3 +a6 b1 b2 b3 −a5 b1 b23 −b2 b23 b5 +b22 b3 b6 +b1 b2 b3 b8 − 2b1 b22 b9 −b33 c5 +b2 b23 c6 +b1 b23 c8 −2b1 b2 b3 c9 −b3 b6 ω 2 +2b1 b9 ω 2 ))+a3 (a22 (a9 b31 −a6 b21 b3 + a4 b1 b23 +2b2 b23 b4 −2b1 b23 b5 +2b21 b3 b8 −2b21 b2 b9 −b33 c4 +b1 b23 c6 −b21 b3 c9 )+2a2 (a6 b1 b22 b3 + a7 b21 b23 − a5 b1 b2 b23 − b22 b23 b5 + b32 b3 b6 + b1 b2 b23 b7 − b1 b32 b9 − b2 b33 c5 + b22 b23 c6 + b1 b33 c7 − b1 b22 b3 c9 − b23 b5 ω 2 + b2 b3 b6 ω 2 + b1 b3 b8 ω 2 − b1 b2 b9 ω 2 + b23 c6 ω 2 − b1 b3 c9 ω 2 + a9 b21 (−b22 + ω 2 ))+ω 2 (a9 b1 (b22 +ω 2 )−b3 (2a8 b1 b2 −3a7 b1 b3 +a5 b2 b3 +b23 c7 −b2 b3 c8 +b22 c9 +c9 ω 2 − a6 (b22 +ω 2 )))))))/(2(a23 b1 +a3 (−a1 +b2 )b3 −a2 b23 )ω 3 (a23 b22 (b22 +ω 2 )+2a2 a3 b2 (a3 b1 b2 − b3 (a1 b2 + b22 + ω 2 )) + a22 (a23 b21 − 2a3 b1 (a1 + b2 )b3 + b23 (a21 + 2a1 b2 + b22 + ω 2 )))) and F = (−a41 b3 (a9 b22 +b3 (−a8 b2 +a7 b3 ))+a43 b1 (b22 c4 −b1 b2 c5 +b21 c7 )+a31 (a2 b3 (2a9 b1 b2 + b3 (−a8 b1 − a6 b2 + a5 b3 − b3 b7 + b2 b8 ) − b22 b9 ) + a3 (a9 b1 b22 − b3 (2a8 b1 b2 − a6 b22 − 3a7 b1 b3 + b3 (a5 b2 + b3 c7 − b2 c8 ) + b22 c9 ))) + a22 b3 (−a2 (a9 b31 − a6 b21 b3 + a4 b1 b23 + b2 b23 b4 − b1 b2 b3 b6 + b21 b2 b9 + b33 c4 − b1 b23 c6 + b21 b3 c9 ) − (a9 b21 + b3 (−a6 b1 + a4 b3 ))ω 2 ) + a33 (a7 b31 b2 − a5 b21 b22 + a4 b1 b32 + a2 b1 b22 b4 + b42 b4 − a2 b21 b2 b5 − b1 b32 b5 + a2 b31 b7 + b21 b22 b7 − 2a2 b1 b2 b3 c4 + b32 b3 c4 + a2 b21 b3 c5 − b1 b22 b3 c5 + a2 b21 b2 c6 + b21 b2 b3 c7 − a2 b31 c8 + (b22 b4 − b1 b2 b5 + b21 b7 )ω 2 ) − a23 (a22 b1 (b2 (2b3 b4 − b1 b6 ) − b23 c4 + b1 b3 (−b5 + c6 ) + b21 (b8 − c9 )) + (a7 b21 +b2 (−a5 b1 +a4 b2 ))b3 ω 2 +a2 (a8 b31 b2 −a6 b21 b22 +a7 b31 b3 −2a5 b21 b2 b3 +3a4 b1 b22 b3 + 3b32 b3 b4 − 2b1 b22 b3 b5 − b1 b32 b6 + b21 b2 b3 b7 + b21 b22 b8 + 3b22 b23 c4 − 2b1 b2 b23 c5 − b1 b22 b3 c6 + b21 b23 c7 + b21 b2 b3 c8 + (2b2 b3 b4 − b1 b3 b5 − b1 b2 b6 + b21 b8 )ω 2 )) + a2 a3 (a22 b1 (b23 b4 − b1 b3 b6 + b21 b9 ) + b3 (a8 b21 − a6 b1 b2 − a5 b1 b3 + 2a4 b2 b3 )ω 2 + a2 (a9 b31 b2 + a8 b31 b3 − 2a6 b21 b2 b3 − a5 b21 b23 + 3a4 b1 b2 b23 + 3b22 b23 b4 − b1 b2 b23 b5 − 2b1 b22 b3 b6 + b21 b2 b3 b8 + b21 b22 b9 + 3b2 b33 c4 − b1 b33 c5 − 2b1 b2 b23 c6 + b21 b23 c8 + b21 b2 b3 c9 + (b23 b4 − b1 b3 b6 + b21 b9 )ω 2 )) + a1 (a33 (a7 b31 − a5 b21 b2 +a4 b1 b22 −b22 b3 c4 +2b1 b2 b3 c5 −b1 b22 c6 −3b21 b3 c7 +b21 b2 c8 )+a2 b3 (−a2 (−2a9 b21 b2 + a8 b21 b3 + a6 b1 b2 b3 − a5 b1 b23 + a2 b23 b4 − b2 b23 b5 − a2 b1 b3 b6 + b22 b3 b6 + b1 b2 b3 b8 + a2 b21 b9 − 2b1 b22 b9 −b33 c5 +b2 b23 c6 +b1 b23 c8 −2b1 b2 b3 c9 )+(2a9 b1 b2 +b3 (−a8 b1 −a6 b2 +a5 b3 ))ω 2 )+ a23 (b2 (a8 b21 b2 − a6 b1 b22 − 2a7 b21 b3 + a5 b1 b2 b3 + b22 b3 b5 − b32 b6 − 2b1 b2 b3 b7 + b1 b22 b8 + b2 b23 c5 − b22 b3 c6 − 2b1 b23 c7 + b1 b2 b3 c8 ) + a2 (−a8 b31 + a6 b21 b2 + a5 b21 b3 − 2a4 b1 b2 b3 − b22 b3 b4 +2b1 b2 b3 b5 −b1 b22 b6 −3b21 b3 b7 +b21 b2 b8 +2b2 b23 c4 −2b1 b23 c5 +2b21 b3 c8 −2b21 b2 c9 )+ (b2 b3 b5 − b22 b6 − 2b1 b3 b7 + b1 b2 b8 )ω 2 ) + a3 (a22 (a9 b31 − a6 b21 b3 + a4 b1 b23 + 2b2 b23 b4 − 2b1 b23 b5 + 2b21 b3 b8 − 2b21 b2 b9 − b33 c4 + b1 b23 c6 − b21 b3 c9 ) − 2a2 (a9 b21 b22 − a6 b1 b22 b3 − a7 b21 b23 + a5 b1 b2 b23 + b22 b23 b5 − b32 b3 b6 − b1 b2 b23 b7 + b1 b32 b9 + b2 b33 c5 − b22 b23 c6 − b1 b33 c7 + LIMIT CYCLES BIFURCATING OF KOLMOGOROV SYSTEMS 13 b1 b22 b3 c9 ) + b3 (−a8 b1 b2 + a6 b22 + 2a7 b1 b3 − a5 b2 b3 )ω 2 + a2 (−b23 b5 + b2 b3 b6 + b1 b3 b8 − 2b1 b2 b9 )ω 2 )) + a21 (a23 (a8 b21 b2 − a6 b1 b22 + b3 (−3a7 b21 + 2a5 b1 b2 − a4 b22 − b2 b3 c5 + b22 c6 + 3b1 b3 c7 − 2b1 b2 c8 ) + b1 b22 c9 ) − a2 b3 (a9 b1 b22 − a8 b1 b2 b3 + a7 b1 b23 + b2 b23 b7 − b22 b3 b8 + b32 b9 + a2 (a9 b21 − a6 b1 b3 + b3 (a4 b3 − b3 b5 + b2 b6 + b1 b8 ) − 2b1 b2 b9 ) + b33 c7 − b2 b23 c8 + b22 b3 c9 ) − b3 (a9 b22 + b3 (−a8 b2 + a7 b3 ))ω 2 + a3 (a9 b1 b32 + a2 (−2a9 b21 b2 + 2a8 b21 b3 − 2a5 b1 b23 + 2a4 b2 b23 − b2 b23 b5 + b22 b3 b6 + 3b1 b23 b7 − 2b1 b2 b3 b8 + b1 b22 b9 + b33 c5 − b2 b23 c6 − b1 b23 c8 + 2b1 b2 b3 c9 ) + b2 (−a8 b1 b2 b3 + a7 b1 b23 + b2 b23 b7 − b22 b3 b8 + b32 b9 + b33 c7 − b2 b23 c8 + b22 b3 c9 ) + (b23 b7 − b2 b3 b8 + b22 b9 )ω 2 )))/((a23 b1 + a3 (−a1 + b2 )b3 − a2 b23 )ω 5 ) Acknowledgements The first author is partially supported by the Ministerio de Economı́a, Industria y Competitividad, Agencia Estatal de Investigación grants MTM2016 - 77278 - P (FEDER) and MDM - 2014 - 0445, the Agència de Gestió d’Ajuts Universitaris i de Recerca grant 2017SGR1617, and the H2020 European Research Council grant MSCA-RISE-2017-777911. The second author is supported by CONICYT PCHA / Postdoctorado en el extranjero Becas Chile / 2018 - 74190062. The third author is partially supported by FCT/Portugal through UID/ MAT/ 04459/2013. References 1. A. Buică and J. Llibre Averaging metods for finding periodic orbits via Brouwer degree, Bull. Sci. Math. 28 (2004). 2. F.H. Busse Transition to turbulence via the statistical limit cycle rout, Syneretics, Springer- Verlag, Berlin, 1978. 3. R.H. Hering, Oscillations in Lotka-Volterra systems of chemical reactions, J. Math. Chem. 5 (1990), 197–202. 4. A. Kolmogorov Sulla teoria di Volterra della lotta per l’esistenza, G. Ist. Ital. Degli Attuari 7 (1936), 74–80. 5. N.G. Lloyd, Degree theory, Cambridge Trends in Mathematics, 73, Cambridge University Press, Cambridge 1978. 6. A.J. Lotka, Elements of physical biology, Science Progress in the Twentieth Century (1919- 1933) 82 (1926), 341–343. 7. G. Laval and R. Pellat Plasma physics. Proceedings of Summer School of Theoretical Physics, Gordon and Breach, New York, 1975. 8. R. M. May Stability and complexity in model ecosystems, Princeton, New Jersey, 1974. 9. R.M. May and W.J. Leonard Nonlinear aspects of competition between three species, SIAM J. Appl. Math. 29 (1975), 243–253. 10. M. Moniri ans S. Moniri, Limit cycles and their period detection via numeric and symbolic hybrid computations, Commun. Nonlinear Sci. Numer. Simul. 83 (2020), 105107, 18 pp. 11. A.V. Razgulin and S.V. Sazonova, Hopf bifurcation in diffusive model of nonlinear optical system with matrix Fourier filtering, Commun. Nonlinear Sci. Numer. Simul. 77 (2019), 288– 304. 12. S. Solomon and P. Richmond, Stable power laws in variable economies; Lotka-Volterra implies Pareto-Zipf, The European Physical Journal B–Condensed Matter and Complex Systems 27 (2002), 257–261. 13. V. Volterra, Variazioni e fluttuazioni del numero d’individui in specie animali conviventi, Memoire della R. Accademia Nazionale dei Lincei II (1926), 558–560. 14 J. LLIBRE, Y.P. MARTÍNEZ AND C. VALLS 1 Departament de Matemàtiques, Facultat de Ciències Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Catalonia, Spain Email address: jllibre@mat.uab.cat 2 Centre de Recerca Matemática, 08193 Bellaterra, Barcelona, Catalonia, Spain Email address: yohanna.martinez@uab.cat 3 Departamento de Matemática, Facultad de Ciencias, Universidad del Bı́o-Bı́o, Casilla 5-C, Concepción, Chile Email address: ymartinez@ubiobio.cl 4 Center for Mathematical Analysis, Geometry and Dynamical Systems, Departa- mento de Matemática, Instituto Superior Técnico, Universidade de Lisboa, 1049-001 Lisboa, Portugal Email address: cvalls@math.tecnico.ulisboa.pt