This is a preprint of: “On the zero-Hopf bifurcatIon of the Lotka-Volterra systems in R3 ”, Maoan Han, Jaume Llibre, Yun Tian, Mathematics, vol. 8, 1137:1–14, 2020. DOI: [10.3390/math8071137] 1 ON THE ZERO–HOPF BIFURCATiON OF THE 2 LOTKA–VOLTERRA SYSTEMS IN R3 3 Maoan Han∗ , Jaume Llibre∗∗1 and Yun Tian∗ ∗ 4 Department of Mathematics, Shanghai Normal University, 100 Guilin Roda, 5 Xuhui District, 200234, Shanghai, China. ∗∗ 6 Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bel- 7 laterra, Barcelona, Catalonia, Spain. Abstract. Here we study the Lotka-Volterra systems in R3 , i.e. the differential systems of the form ∑3 dxi ( ) = xi ri − aij xj , i = 1, 2, 3. dt j=1 8 It is known that some of these differential systems can have at least four periodic 9 orbits bifurcating from one of their equilibrium points. Here we prove that there are 10 some of these differential systems exhibiting at least six periodic orbits bifurcating 11 from one of their equilibrium points. The tool for proving this result is the averaging 12 theory of third order. 13 2100 Mathematics Subject Classification. Primary 34C07, 34C08, 37G15. 14 Key words and phrases. Lotka–Volterra polynomial differential systems, periodic 15 orbit, Hopf bifurcation, averaging theory. 1Jaume Llibre is the corresponding author 1 http://www.gsd.uab.cat 2 1 1. Introduction and statement of results 2 An equilibrium point of a 3-dimensional autonomous differential system having a 3 pair of purely imaginary eigenvalues and a zero eigenvalue is a zero-Hopf equilibrium. 4 A 2-parameter unfolding of a 3-dimensional autonomous differential system with 5 a zero-Hopf equilibrium is a zero-Hopf bifurcation. More precisely, when the two 6 parameters of the unfolding are zero we have an isolated zero-Hopf equilibrium, 7 and the dynamics of the unfolding is complex and sometimes chaotic in a small 8 neighborhood of this isolated equilibrium when we vary the two parameters in a 9 small neighborhood of the origin, see for more details [4, 7, 8, 17, 26] and references 10 quoted there. 11 A Lotka-Volterra system in R3 with coordinates (x1 , x2 , x3 ) is a quadratic poly- 12 nomial differential system of the form ∑3 dxi ( ) (1) = xi ri − aij xj , i = 1, 2, 3, dt j=1 13 where the dot denotes derivative with respect to the independent variable t, usually 14 called the time, and the ri ’s and the aj ’s are parameters. 15 Many natural phenomena can be modeled by the Lotka–Volterra systems, start- 16 ing in biology with the time evolution of conflicting species that now continuing 17 being studied intensively see [9, 10, 11, 12, 13, 14, 15, 22, 25, 27, 28, 31], later on 18 problems of plasma physics [18], or problems in hydrodynamics [3], ... 19 It is known that Lotka-Volterra systems can exhibit zero-Hopf equlibria, see for 20 instance [20]. Then a natural question is if we perturbed a Lotka-Volterra system (1) 21 having a zero-Hopf equilibrium point inside the class of all Lotka-Volterra systems 22 how many periodic orbits can bifurcate from such an equilibrium? 23 Note that the unfolding of Lotka-Volterra system (1) with a zero-Hopf equilib- 24 rium needs at least a 3-parameter family. Arnold [1] in 1973 proposed to investigate 25 bifurcations of 3-parameter families with a zero–Hopf equilibrium. 26 As far as we know the number of periodic orbits which can bifurcate from a zero- 27 Hopf equilibrium point when this is perturbed inside the class of all Lotka-Volterra 28 systems only has been studied partially in the paper [20] using averaging theory of 29 second order. There the authors provided explicit conditions for the existence of 30 one or two periodic orbits bifurcating from one of these equilibria. 31 Here we shall use the averaging theory of third order for studying the num- 32 ber of periodic orbits which can bifurcate from a zero-Hopf equilibrium point of 33 a Lotka-Volterra system (1). Previous results in this direction are the following. 34 First we say that an equilibrium point of a 3-dimensional autonomous differential 35 system having a pair of purely imaginary eigenvalues and a non-zero eigenvalue is 36 a Hopf equilibrium. The bifurcation of periodic orbits in a Hopf equilibrium of a 37 Lotka-Volterra system (1) have been studied by many authors. Thus in the papers 38 [16, 23, 30] the authors proved that two periodic orbits can bifurcate from a Hopf 39 equilibrium of system (1). While in [5, 6, 24] it is shown that three periodic orbits 40 can bifurcate from a Hopf equilibrium. Recently in [29] it is proved that four peri- 41 odic orbits can bifurcate from a Hopf equilibrium of system (1). All these previous 3 1 results on the number of periodic orbits bifurcating from a Hopf equilibrium are 2 when system (1) has all its coefficients aij and ri positive, and under this assump- 3 tion in [5] it is conjectured that at least five periodic orbits can bifurcate from a 4 such Hopf equilibrium, but this conjecture remains open. 5 In short, until now it is known that there are Lotka-Volterra systems (1) having 6 at least four periodic orbits bifurcating from one of their equilibrium points. Our 7 main result is the following one. 8 Theorem 1. There are Lotka-Volterra systems (1) having at least six periodic 9 orbits bifurcating from a zero–Hopf equilibrium. 10 We remark that those Lotka-Volterra systems (1) exhibiting a Hopf bifurcation 11 with at least six periodic orbits do not have all the coefficients aij and ri positive. 12 The proof of Theorem 1 is given in the next section. 13 2. Proof of Theorem 1 14 If system (1) has a zero-Hopf equilibrium (a, b, c) with non-zero components 15 without loss of generality we can consider this equilibrium at the point (1, 1, 1) 16 doing the rescaling (x, y, z) → (x/a, y/b, z/c). Then every Lotka-Volterra system 17 (1) having the equilibrium (1, 1, 1) can be written as ( ) ẋ = x a11 (x − 1) + a12 (y − 1) + a13 (z − 1) , ( ) (2) ẏ = y a21 (x − 1) + a22 (y − 1) + a23 (z − 1) , ( ) ż = z a31 (x − 1) + a32 (y − 1) + a33 (z − 1) , 18 where now we denote the coordinates of R3 by (x, y, z). Since we shall use the 19 averaging theory of third order for studying the periodic orbits of this system we 20 take the coefficients aij as follows aij = aij0 + εaij1 + ε2 aij2 + ε3 aij3 , 21 with i and j varying in {1, 2, 3}, being ε a small parameter. Note that in the 22 differential system (2) there are 37 parameters. This big number of parameters 23 produce that the computations for studying the number of periodic orbits which can 24 bifurcate from the equilibrium (1, 1, 1) are tedious and huge. All the computations 25 of this paper has been done with the help of the algebraic manipulator mathematica. 26 First we translate the equilibrium (1, 1, 1) to the origin of coordinates and system 27 (2) becomes ( ẋ = (1 + x) a110 x + a120 y + a130 z + ε(a111 x + a121 y + a131) z)+ ε2 (a112 x + a122 y + a132 z) + ε3 (a113 x + a123 y + a133 z) , ( ẏ = (1 + y) a210 x + a220 y + a230 z + ε(a211 x + a221 y + a231)z)+ (3) ε2 (a212 x + a222 y + a232 z) + ε3 (a213 x + a223 y + a233 z) , ( ż = (1 + z) a310 x + a320 y + a330 z + ε(a311 x + a321 y + a331)z)+ ε2 (a312 x + a322 y + a332 z) + ε3 (a313 x + a323 y + a333 z) . 28 Choosing the conditions (4) a110 = a120 = a130 = a210 = 0, a320 = −(a2220 + ω 2 )/a230 and a330 = −a220 , 4 1 with a230 ω ̸= 0 it is easy to check that the linear part of system (3) at the origin 2 has eigenvalues 0 and ±ωi. So the origin of system (3) is a zero-Hopf equilibrium, 3 and consequently system (2) has a zero-Hopf equilibrium at the point (1, 1, 1). We 4 remark that there are other conditions which also provide that the point (1, 1, 1) 5 be a zero-Hopf equilibrium. 6 In what follows we shall study the periodic orbits bifurcating from the zero-Hopf 7 equilibrium (0, 0, 0) of system (3) under conditions (4). 8 As we shall see the amount of computations for studying this Hopf-bifurcation 9 are huge due to the big number of parameters in system (3). 10 In order to study the periodic orbits bifurcating from the zero-Hopf equilibrium 11 at the origin of the differential system (3) using the averaging theory of third 12 order (see the appendix), we need to introduce a small parameter and take a new 13 independent variable in which the differential system be periodic. 14 The small parameter for the averaging theory will be the parameter ε, and we 15 do the rescaling (x, y, z) = (εX, εY, εZ). Then system (3) in the new variables 16 (X, Y, Z) writes Ẋ = ε(a111 X + a121 Y + a131 Z) + ε2 (a112 X + a111 X 2 + a122 Y + a121 XY + a132 Z + a131 XZ) + ε3 (a113 X + a112 X 2 + a123 Y + a122 XY + a133 Z + a132 XZ) + O(ε4 ), Ẏ = a220 Y + a230 Z + ε(a211 X + a221 Y + a220 Y 2 + a231 Z + a230 Y Z)+ ε2 (a212 X + a222 Y + a211 XY + a221 Y 2 + a232 Z + a231 Y Z)+ (5) ε3 (a213 X + a223 Y + a212 XY + a222 Y 2 + a233 Z + a232 Y Z) + O(ε4 ), Ż = (a230 a310 X − a2220 Y − a220 a230 Z − Y ω 2 )/a230 + ε(a230 a311 X+ a230 a321 Y + a230 a331 Z + a230 a310 XZ − a2220 Y Z − a220 a230 Z 2 − Y Zω 2 )/a230 + ε2 (a312 X + a322 Y + a332 Z + a311 XZ + a321 Y Z+ a331 Z 2 ) + ε3 (a313 X + a323 Y + a333 Z + a312 XZ + a322 Y Z+ a332 Z 2 ) + O(ε4 ). In order to simplify the computations of the averaging theory we shall write the linear part of the differential system (5) into its real Jordan normal form doing the linear change of variables (X, Y, Z) → (u, v, w) given by X = w, a230 a310 w a230 ωv − a220 a230 u Y = + , ω2 a2220 + ω 2 Z = −a220 a230 a310 w + a230 ω 2 u. 5 1 Now the differential system (5) in the new variables (u, v, w) becomes ε ( u̇ = −ωv + (a131 a3220 a310 − a121 a2220 a230 a310 + ω 4 (a2220 + ω 2 ) a131 a220 a310 ω 2 − a220 a230 a321 ω 2 + a2220 a331 ω 2(+ a331 ω 4 )ω 2 u+ a230 (a121 a220 a310 + a321 ω 2 )ω 3 v − (a2220 + ω 2 ) (a131 a2220 a2310 − a121 a220 a230 a2310 − a111 a220 a310 ω 2 − a230 a310 a321 ω 2))+ a220 a310 a331 ω 2 − a311 ω 4 )w − ω 5 uv + a220 a310 ω 3 vw + O(ε2 ), ε ( v̇ = ωu + (−a3220 a221 a230 + a4220 a231 − a230 ω 3 (a2220 + ω 2 ) a131 a2220 a230 a310 + a121 a220 a2230 a310 − a2220 a2230 a321 + a3220 a230 a331 − a220 a221 a230 ω 2 + 2a2220 a231 ω 2 − a131 a230 a310 ω 2 + a220 a230 a331 ω 2 + a231 ω 4 )ω 2 u + a230 (a2220 a221 − a121 a230 a310 + (6) a220 a230 a321 + a221 ω 2 )ω 3 v − (a2220 + ω 2 )(−a2220 a221 a230 a310 + a3220 a231 a310 − a131 a220 a230 a2310 + a121 a2230 a2310 − a220 a2230 a310 a321 + a2220 a230 a310 a331 − a211 a2220 ω 2 + a111 a230 a310 ω 2 − a221 a230 a310 ω 2 + a220 a231 a310 ω 2 − a220 a230 a311 ω 2 − a211 ω 4 )w − a220 a2230 ω 4 u2 − a230 (a3220 + a2220 a230 + a220 ω 2 − a230 ω 2 )ω 3 uv + a2230 a310 (a2220 + ω 2 )ω)2 uw+ a220 a2230 v 2 ω 4 + a220 a230 (a220 + a230 )a310 (a2220 + ω 2 )ωvw + O(ε2 ), ε ( ẇ = (a131 a2220 − a121 a220 a230 + a131 ω 2 )ω 2 u+ ω 2 (a2220 + ω 2 ) ) a121 a230 ω 3 v − (a2220 + ω 2 )(a131 a220 a310 − a121 a230 a310 − a111 ω 2 )w +O(ε2 ). 2 In the computations of the previous differential system we have obtained the ex- 3 pressions of u̇, v̇ and ẇ until terms of O(ε4 ), but here we only present them until 4 terms of order O(ε2 ), otherwise the expression of system (6) would need several 5 pages. Using an algebraic manipulator as mathematica or mapple it is relatively 6 easy to repeat our computations. 7 Now we write the differential system (6) in cylindrical coordinates (r, θ, w) where 8 u = r cos θ and v = r sin θ, and taking θ as the new independent variable of the 9 differential system defined we get the new differential system r′ = εF11 (θ, r, w) + ε2 F21 (θ, r, w) + ε3 F31 (θ, r, w) + O(ε4 ), (7) w′ = εF12 (θ, r, w) + ε2 F22 (θ, r, w) + ε3 F32 (θ, r, w) + O(ε4 ), defined in in r > 0, where the prime denotes derivative with respect to the variable θ. Here we only provide the explicit expressions of F11 = F11 (θ, r, w) and F12 = F12 (θ, r, w) which are the shorter ones, but our next computations will use the 6 expressions of F21 , F22 , F31 and F32 . Thus we have 1 ( F11 = 5 2 2 (a230 (a131 a3220 a310 − a121 a2220 a230 a310 + a230 ω (a220 + ω ) a131 a220 a310 ω 2 − a220 a230 a321 ω 2 + a2220 a331 ω 2 + a331 ω 4 ) cos2 θ+ (a4220 a231 − a3220 a221 a230 − a131 a2220 a230 a310 + 2a121 a220 a2230 a310 − a2220 a2230 a321 + a3220 a230 a331 − a220 a221 a230 ω 2 + 2a2220 a231 ω 2 − a131 a230 a310 ω 2 + a2230 a321 ω 2 + a220 a230 a331 ω 2 + a231 ω 4 )ω cos θ sin θ+ a230 (a2220 a221 − a121 a230 a310 + a220 a230 a321 + a221 ω 2 )ω 2 sin2 θ)ω 2 r− (a2220 + ω 2 )(a230 (a131 a2220 a2310 − a121 a220 a230 a2310 − a111 a220 a310 ω 2 − a230 a310 a321 ω 2 + a220 a310 a331 ω 2 − a311 ω 4 ) cos θ + (a3220 a231 a310 − a2220 a221 a230 a310 − a131 a220 a230 a2310 + a121 a2230 a2310 − a220 a2230 a310 a321 + a2220 a230 a310 a331 − a211 a2220 ω 2 + a111 a230 a310 ω 2 − a221 a230 a310 ω 2 + a220 a231 a310 ω 2 − a220 a230 a311 ω 2 − a211 ω 4 )ω sin θ])w − (a230 ω(a2220 + a220 a230 + ω 2 ) cos2 θ sin θ + a230 (a3220 + a2220 a230 + a220 ω 2 − a230 ω 2 ) cos θ sin2 θ − a220 a2230 ω sin3 θ)ω 4 r2 + (a230 (a220 + a230 )a310 ω 3 (a2220) + ω 2 ) cos θ sin θ + a220 a230 (a220 + a230 )a310 ω 2 (a2220 + ω 2 ) sin2 θ)rw , 1 ( F12 = 3 2 2 ((a131 a2220 − a121 a220 a230 + a131 ω 2 ) cos θ+ ω (a220 + ω ) a121 a230 ω)sin θ)ω 2 r − (a2220 + ω 2 )(a131 a220 a310 − a121 a230 a310 − a111 ω 2 )w , We note that the differential system (7) is written in the normal form (11) for applying the averaging theory of third order described in the appendix, where the variables t and x of the appendix are now θ and (r, w) respectively. Computing the averaged function of first order f1 (r, w) = (f11 (r, w), f12 (r, w)) defined in the appendix we get f11 (r, w) = Ar, f12 (r, w) = Bw, where (a131 a220 − a121 a230 )a310 + (a221 + a331 )ω 2 + (a220 + a230 )a310 a220 w A= , 2ω 3 (a121 a230 − a131 a220 )a310 + a111 ω 2 B= . ω3 1 We look for the zeros (r∗ , w∗ ) of f1 (r, w) with r > 0, and since the unique zero of the 2 function f1 (r, w) is (0, 0), or a continuum of zeros if the coefficient A or B is zero, 3 the averaged function of first order does not give any information on the periodic 4 solutions of system (7), see the appendix. Therefore we force that the averaged 5 function of first order be identically zero and we shall use the averaged functions 6 of higher order to obtain information on the periodic solutions of the differential 7 system (7). 8 Since the coefficient of rw in the function f11 (r, w) is (a220 + a230 )a310 a220 we 9 need to consider the following three cases in order that the averaged function of 10 first order be identically zero: 11 Case 1: a220 = −a230 , 12 a331 = (a121 a230 a310 + a131 a230 a310 − a221 ω 2 )/ω 2 , 13 a111 = (−a121 a230 a310 − a131 a230 a310 )/ω 2 . 7 1 Case 2: a310 = 0, a331 = −a221 , a111 = 0. 2 Case 3: a220 = 0, 3 a331 = (a121 a230 a310 − a221 ω 2 )/ω 2 , 4 a111 = −(a121 a230 a310 )/ω 2 . Case 1. Since the averaged function of first order f1 (r, w) is identically zero, we compute the averaged function of second order f2 (r, w) = (f21 (r, w), f22 (r, w)) and we obtain f21 (r, w) = (Cw + D)r, f22 (r, w) = Ew, where C = −(−a121 a2230 a2310 − a131 a2230 a2310 + a121 a230 a310 ω 2 + a131 a230 a310 ω 2 + a221 a230 a310 ω 2 + a230 a231 a310 ω 2 + a211 ω 4 )/(2ω), D = −(a2121 a4230 a2310 + 2a121 a131 a4230 a2310 + a2131 a4230 a2310 − 2a121 a221 a3230 a310 ω 2 −2a131 a221 a3230 a310 ω 2 − a121 a3230 a231 a310 ω 2 − a131 a2303 a231 a310 ω 2 + a121 a230 a310 a321 ω + a131 a230 a310 a321 ω − a121 a211 a230 ω 4 − 3 2 3 2 2 a131 a211 a2230 ω 4 − a131 a221 a230 a310 ω 4 + a122 a2230 a310 ω 4 + a132 a2230 a310 ω 4 −a131 a230 a231 a310 ω 4 + a121 a2230 a311 ω 4 + a131 a2230 a311 ω 4 − a131 a211 ω 6 − a222 a230 ω 6 − a230 a332 ω 6 )/(2a230 ω 7 ), E = (a2121 a4230 a2310 + 2a121 a131 a4230 a2310 + a2131 a4230 a2310 − 2a121 a221 a3230 a310 ω 2 −2a131 a221 a3230 a310 ω 2 − a121 a3230 a231 a310 ω 2 − a131 a3230 a231 a310 ω 2 + a121 a3230 a310 a321 ω 2 + a131 a3230 a310 a321 ω 2 − a121 a211 a2230 ω 4 − a131 a211 a2230 ω 4 − a131 a221 a230 a310 ω 4 + a122 a2230 a310 ω 4 + a132 a2230 a310 ω 4 −a131 a230 a231 a310 ω 4 + a121 a2230 a311 ω 4 + a131 a2230 a311 ω 4 − a131 a211 ω 6 + a112 a230 ω 6 )/(a230 ω 7 ). Again the unique zero of the averaged function of second order f2 (r, w) is the (0, 0) or a continuum of solutions in case that convenient coefficients C, D or E are zero. Therefore the averaging theory of second order does not provide any information on the periodic solutions of the differential system (7). Consequently we impose that the averaged function of second order f2 (r, w) be identically zero, and we obtain that ( a211 = a121 a2230 a2310 + a131 a2230 a2310 − a121 2 ) a4230 a310 ω − a131 a230 a310 ω − 2 2 2 a221 a230 a310 ω − a230 a231 a310 ω /ω , ( 2 2 2 a332 = a121 a230 a310 + 2a121 a131 a230 a310 + a2131 a2230 a310 − a121 a221 a2230 a310 − a131 a221 a2230 a310 − a121 a131 a230 a2310 − a2131 a230 a2310 + a121 a2230 a310 a321 + a131 a2230 a310 a321 + a121 a131 a310 ω 2 + a2131 a310 ω 2 + a122 a230 a)310 ω 2 + a132 a230 a310 ω 2 + a121 a230 a311 ω 2 + a131 a230 a311 ω 2 − a222 ω 4 /ω 4 , ( a112 = − a2121 a2230 a310 − 2a121 a131 a2230 a310 − a2131 a2230 a310 + a121 a221 a2230 a310 +a131 a221 a2230 a310 + a121 a131 a230 a2310 + a2131 a230 a2310 − a121 a2230 a310 a321 −a131 a2230 a310 a321 − a121 a131 a310 ω 2 − a2131 a310 ω 2 − ) a122 a230 a310 ω 2 − 2 2 2 4 a132 a230 a310 ω − a121 a230 a311 ω − a131 a230 a311 ω /ω . 8 We compute the averaged function of third order f3 (r, w) = (f31 (r, w), f32 (r, w)) and we get a0 r4 + a1 r3 + a2 r2 w + a3 r2 + a4 rw + a5 w2 + a6 r + a7 w f31 (r, w) = , 384a230 (a2230 + ω 2 )ω 13 r b0 r3 + b1 r2 w + b2 r2 + b3 rw + b4 r + b5 w f32 (r, w) = . 24a230 (a2230 + ω 2 )2 ω 9 1 We do not provide the explicit expressions of the coefficients aj and bj because we 2 shall need approximately twenty pages for writing them. 3 Now we shall study the zeros of the function f3 (r, w). Since the variable w 4 appears linearly in the equation f32 (r, w) = 0, we isolate it and we get w = W (r). 5 Substituting w = W (r) into the equation f31 (r, w) = 0, we obtain an equation in 6 the variable r of the form n(r) c2 r2 + c3 r3 + c4 r4 + c5 r5 + c6 r6 + c7 r7 + c8 r8 (8) = = 0. d(r) (d0 + d1 r + d2 r2 )2 7 The coefficients cj and dj are polynomials in some of the coefficients of the differen- 8 tial system (2), more precisely in the coefficients a113 , a121 , a122 , a123 , a131 , a132 , a133 , 9 a212 , a221 , a222 , a223 , a230 , a231 , a232 , a310 , a311 , a312 , a321 , a322 , a333 , ω. We have com- 10 puted the rank of the Jacobian matrix of the function (c2 , c3 , c4 , c5 , c6 , c7 , c8 ) with 11 respect to the 21 previous coefficients, it is the rank of a 7 × 23 matrix, and we get 12 that this rank is 7. Therefore the seven coefficients of the polynomial n(r) are inde- 13 pendent, and consequently we can choose them in such a way that the polynomial 14 n(r) has six positive real roots. Moreover, we also can choose those coefficients in 15 such a way that the resultant of the polynomials n(r) and d(r) is not zero, and 16 consequently both polynomials do not have a common root. So equation (8) can 17 have six positive solutions, rj∗ for j = 1, 2, 3, 4, 5, 6. 18 In short, we have that (rj∗ , W (rj∗ )) for j = 1, 2, 3, 4, 5, 6 are six zeros of the third 19 averaged function f3 (r, w). These zeros can be chosen simple, i.e. the Jacobian 20 of the function f3 (r, w) evaluated in such zeros is not zero. Consequently by the 21 averaging theory (see the appendix) the differential system (7) has six periodic 22 solutions (rj (θ, ε), wj (θ, ε)) such that (rj (0, ε), wj (0, ε)) → (rj∗ , W (rj∗ )) when ε → 0. Going back to the differential system (6) we obtain for this system six periodic solutions (uj (t, ε), vj (t, ε), wj (t, ε)) such that (uj (0, ε), vj (0, ε), wj (0, ε)) → (rj∗ , 0, W (rj∗ )), when ε → 0. These periodic solutions provide six periodic solutions (Xj (t, ε), Yj (t, ε), Zj (t, ε)) for the differential system (5) such that Xj (0, ε) → W (rj∗ ), a230 a310 W (rj∗ ) a220 a230 rj∗ Yj (0, ε) → − 2 , ω2 a220 + ω 2 Zj (0, ε) → a230 ω 2 rj∗ − a220 a230 a310 W (rj∗ ), 9 1 when ε → 0. Finally going back to the differential system (2) we obtain six periodic 2 solutions (xj (t, ε), yj (t, ε), zj (t, ε)) such that xj (0, ε) = 1 + εW (rj∗ ) + O(ε2 ), ( a230 a310 W (rj∗ a220 a230 rj∗ ) (9) yj (0, ε) = 1 + ε 2 − 2 + O(ε2 ), ω a220 + ω 2 ( ) zj (0, ε) = 1 + ε a230 ω 2 rj∗ − a220 a230 a310 W (rj∗ ) + O(ε2 ), 3 when ε → 0. Clearly from (9) these six periodic solutions (xj (t, ε), yj (t, ε), zj (t, ε)) 4 tend to the equilibrium point (1, 1, 1) of the differential system (2) when ε → 0. 5 Hence they bifurcate from that zero-Hopf equilibrium at ε = 0. This completes the 6 proof of Theorem 1. Case 2. Again since the averaged function of first order f1 (r, w) is identically zero, we compute the averaged function of second order f2 (r, w) = (f21 (r, w), f22 (r, w)) and we obtain f21 (r, w) = (Cw + D)r, f22 (r, w) = Ew, where a220 (a211 a2220 + a211 a220 a230 + a220 a230 a311 + a2230 a311 + a211 ω 2 ) C= , 2a230 ω 3 1 ( D= a121 a211 a220 a230 − a131 a211 a2220 − a131 a220 a230 a311 + 2a230 ω 3 ) a121 a2230 a311 − a131 a211 ω 2 − a222 a230 ω 2 − a230 a332 ω 2 , 1 ( E= a121 a211 a220 a230 − a131 a211 a2220 − a131 a220 a230 a311 + a230 ω 3 ) a121 a2230 a311 − a131 a211 ω 2 + a112 a230 ω 2 . 7 As in the previous case the unique zero of the averaged function of second order 8 f2 (r, w) is the (0, 0) or a continuum of solutions in case that convenient coefficients 9 C, D or E are zero. Consequently we impose that the averaged function of second 10 order f2 (r, w) be identically zero, but since the coefficient of rw in the function 11 f21 (r, w) is a product of two factors we have two consider two subcases. Subcase 2.1: a220 = 0. Then in order that the averaged function of second order f2 (r, w) be identically zero we take a121 a2230 a311 − a131 a211 ω 2 − a222 a230 ω 2 a332 = , a230 ω 2 a131 a211 ω 2 − a121 a2230 a311 a112 = . a230 ω 2 12 We compute the averaged function of third order f3 (r, w) = (f31 (r, w), f32 (r, w)) 13 and we get a0 r3 + a1 r2 + a2 rw + a3 w2 + a4 r + a5 w f31 (r, w) = , 384a2230 ω 5 r (10) b0 r3 + b1 r2 + b2 rw + b3 r + b4 w f32 (r, w) = . 24a2230 ω 5 10 1 Here the expressions of the coefficients aj ’s and bj ’s are relatively short, but we do 2 not need them explicitly. We shall study the zeros of the function f3 (r, w). Since the variable w appears linearly in the equation f32 (r, w) = 0, we isolate it and we get w = W (r). Sub- stituting w = W (r) into the equation f31 (r, w) = 0, we obtain an equation in the variable r of the form c2 r2 + c3 r3 + c4 r4 + c5 r5 + c6 r6 = 0. (d0 + d1 r + d2 r2 )2 3 So at most we have four positive solutions for the variable r, and consequently at 4 most four zeros for the averaged function of third order f3 (r, w). In any case less 5 than the six obtained in Case 1. Subcase 2.2: a211 a2220 + a211 a220 a230 + a220 a230 a311 + a2230 a311 + a211 ω 2 = 0. Then in order that the averaged function of second order f2 (r, w) be identically zero we take a211 a2220 + a211 a220 a230 + a211 ω 2 a311 = − , a230 (a220 + a230 ) a121 a211 + a131 a211 + a220 a222 + a222 a230 a332 = − , a220 + a230 a121 a211 + a131 a211 a112 = . a220 + a230 6 We compute the averaged function of third order f3 (r, w) = (f31 (r, w), f32 (r, w)) 7 and we get again the expressions given in (10), of course the coefficients aj ’s and 8 bj ’s are now different. Repeating the arguments of the previous subcase we obtain 9 at most four zeros for the averaged function of third order f3 (r, w). Case 3. Again since the averaged function of first order f1 (r, w) is identically zero, we compute the averaged function of second order f2 (r, w) = (f21 (r, w), f22 (r, w)) and we obtain f21 (r, w) = (Cw + D)r, f22 (r, w) = Ew, where a310 (a121 − a221 )a230 C= − , 2 ω3 1 ( D= − a121 a3230 a310 a321 − a131 a221 a230 a310 ω 2 + a122 a2230 a310 ω 2 + 2a230 ω 5 ) a121 a2230 a311 ω 2 − a131 a211 ω 4 − a222 a230 ω 4 − a230 a332 ω 4 , 1 ( E= a121 a3230 a310 a321 − a131 a221 a230 a310 ω 2 + a122 a2230 a310 ω 2 + a230 ω 5 ) a121 a2230 a311 ω 2 − a131 a211 ω 4 + a112 a230 ω 4 . 10 As in the previous case the unique zero of the averaged function of second order 11 f2 (r, w) is the (0, 0) or a continuum of solutions in case that convenient coefficients 12 C, D or E are zero. Consequently we impose that the averaged function of second 13 order f2 (r, w) be identically zero, but since the coefficient of rw in the function 14 f21 (r, w) is a product of two factors which can be zero, namely a310 (a121 − a221 ), 15 we have two consider two subcases. 11 Subcase 3.1: a310 = 0. Then in order that the averaged function of second order f2 (r, w) be identically zero we take a121 a2230 a311 − a131 a211 ω 2 − a222 a230 ω 2 a332 = , a230 ω 2 −a121 a2230 a311 + a131 a211 ω 2 a112 = . a230 ω 2 1 We compute the averaged function of third order f3 (r, w) = (f31 (r, w), f32 (r, w)) 2 and we get again the expression given in (10), consequently at most four solutions. Subcase 3.2: a221 = a121 . Then in order that the averaged function of second order f2 (r, w) be identically zero we take 1 ( a332 = a121 a3230 a310 a321 − a121 a131 a230 a310 ω 2 + a122 a2230 a310 ω 2 + a230 ω 2 ) a121 a2230 a311 ω 2 − a131 a211 ω 4 − a222 a230 ω 4 , 1 ( a112 = a121 a131 a230 a310 ω 2 − a121 a3230 a310 a321 − a122 a2230 a310 ω 2 − a230 ω 4 ) a121 a2230 a311 ω 2 + a131 a211 ω 4 . We compute the averaged function of third order f3 (r, w) = (f31 (r, w), f32 (r, w)) and we get a0 r4 w + a1 r4 + a2 r2 w2 + a3 r3 + a4 r2 w + a5 w3 + a6 r2 + a7 rw + a8 w2 f31 (r, w) = , 384a2230 ω 9 r b0 r3 + b1 r2 w + b2 r2 + b3 rw + b4 w2 + b5 r + b6 w f32 (r, w) = . 24a2230 ω 7 3 Here the explicit expressions of the coefficients aj ’s and bj ’s only should need ap- 4 proximately three pages for writing them. But unfortunately in this case we do 5 not know how to control the zeros (r∗ , w∗ ) of the function f3 (r, w) with r∗ > 0. 6 We think that in this subcase it is possible that more than six simple zeros can be 7 obtained, but for the moment this is an open problem. 8 Appendix: The averaging theory of first, second and third order 9 The averaging theory of third order for studying periodic orbits was developed 10 [2] and in [19] at any order. It can be summarized as follows. 11 Consider the differential system (11) ẋ = εF1 (t, x) + ε2 F2 (t, x) + ε3 F3 (t, x) + ε4 R(t, x, ε), 12 where F1 , F2 , F3 : R × D → R, R : R × D × (−εf , εf ) → R are continuous functions, 13 T –periodic in the first variable, and D is an open subset of Rn . Assume that the 14 following hypotheses (i) and (ii) hold. 15 (i) F1 (t, ·) ∈ C 2 (D), F2 (t, ·) ∈ C 1 (D) for all t ∈ R, F1 , F2 , F3 , R, Dx2 F1 , Dx F2 16 are locally Lipschitz with respect to x, and R is twice differentiable with 17 respect to ε. 12 1 We define Fk0 : D → R for k = 1, 2, 3 as ∫ 1 T f1 (x) = F1 (s, x)ds, T 0 ∫ 1 T f2 (x) = [Dx F1 (s, x) · y1 (s, x) + F2 (s, x)] ds, T 0 ∫ 1 T [1 ∂ 2 F1 1 ∂F1 f3 (x) = y1 (s, x)T (s, x)y1 (s, x) + (s, x)y2 (s, x) T 0 2 ∂x2 ] 2 ∂x ∂F2 + (s, x)y1 (s, x) + F3 (s, x) ds, ∂x 2 where ∫ s y1 (s, x) = F1 (t, x)dt, ∫0 s [ ∫ t ] ∂F1 y2 (s, x) = (t, x) F1 (r, x)dr + F2 (t, x) dt. 0 ∂x 0 3 (ii) For an open and bounded set V ⊂ D and for each ε ∈ (−εf , εf ) \ {0}, there 4 exists a ∈ V such that f1 (a) + εf2 (a) + ε2 f3 (a) = 0 and dB (f1 + εf2 + 5 ε2 f3 , V, aε ) ̸= 0 (i.e. the Brouwer degree of the function f1 + εf2 + ε2 f3 at 6 the point a is not zero). 7 Then for |ε| > 0 sufficiently small there exists a T –periodic solution x(t, ε) of system 8 (11) such that x(0, ε) → a when ε → 0. 9 A sufficient condition in order that dB (f1 + εf2 + ε2 f3 , V, aε ) ̸= 0 is that the 10 Jacobian of the function f1 + εf2 + ε2 f3 at a is not zero, see for details [21]. 11 The averaging theory of first order takes place when f1 is not identically zero. 12 Therefore the zeros of f1 + εf2 + ε2 f3 are mainly the zeros of f1 for ε sufficiently 13 small. 14 The averaging theory of second order takes place when f1 is identically zero and 15 f2 is not identically zero. Then the zeros of f1 + εf2 + ε2 f3 are mainly the zeros of 16 f2 for ε sufficiently small. 17 Finally the averaging theory of third order takes place when f1 and f2 are iden- 18 tically zero and f3 is not identically zero. Therefore the zeros of f1 + εf2 + ε2 f3 are 19 mainly the zeros of f3 for ε sufficiently small. 20 Acknowledgments 21 The first and the third authors are partially supported by National Natural 22 Sciences Foundation 11431008, 11771296 and 11501370. 23 The second author is partially supported by the MINECO grants MTM2013- 24 40998-P and MTM2016-77278-P (FEDER) and the AGAUR grant 2014-SGR568. 25 References 26 [1] V.I. Arnold, “Arnold’s problems”, Springer, PHASIS, 2000. 27 [2] A. Buica and J. Llibre, Averaging methods for finding periodic orbits via Brouwer degree, 28 Bull. Sci. Math. 128 (2004), 7–22. 13 1 [3] F.H. Busse “Transition to turbulence via the statistical limit cycle route”, Synergetics, 2 Springer–Verlag, Berlin, p. 39, 1978. 3 [4] A.R. Champneys and V. Kirk, The entwined wiggling of homoclinic curves emerging from 4 saddle-node/Hopf instabilities. Phys. D 195 (2004), 77–105. 5 [5] M. Gyllenberg and P. Yan, On the number of limit cycles for the three-dimensional Lotka- 6 Volterra systems, Discrete Contin. Dyn. Syst. Ser. S 11 (2009), 347–352. 7 [6] M. Gyllenberg, P. Yan and Y. Wang, A 3D competitive LotkaVolterra system with three limit 8 cycles: a falsification of a conjecture by Hofbauer and So, Appl. Math. Lett. 19 (2006), 1–7. 9 [7] J. Guckenheimer, “On a codimension two bifurcation”, Lecture Notes in Math. 898 (1980), 10 99–142. 11 [8] J. Guckenheimer and P. Holmes, “Nonlinear Oscillations, Dynamical Systems and Bifurca- 12 tions of Vector Fields”, Springer, 1983. 13 [9] M.W. Hirsch, Systems of differential equations that are competitive or cooperative. I: Limit 14 sets, SIAM J. Math. Anal. 13 (1982), 167–169. 15 [10] M.W. Hirsch, Systems of differential equations that are competitive or cooperative. II: Con- 16 vergence almost everywhere, SIAM J. Math. Anal. 16 (1985), 423–439. 17 [11] M.W. Hirsch, Stability and convergence in stronly monotone dynamical systems, J. für die 18 reine und angewandte Mathematik 383 (1988), 1–53. 19 [12] M.W. Hirsch, Systems of differential equations that are competitive or cooperative. III: Com- 20 peting species, Nonlinearity 1 (1988), 117–124. 21 [13] M.W. Hirsch, Systems of differential equations that are competitive or cooperative. IV: Struc- 22 tural stability in 3-dimensional systems, SIAM J. Math. Anal. 21 (1990), 1225–1234. 23 [14] M.W. Hirsch, Systems of differential equations that are competitive or cooperative. V: Con- 24 vergence in 3-dimensional systems, J. Differential Equations 80 (1989), 94–106. 25 [15] M.W. Hirsch, Systems of differential equations that are competitive or cooperative. VI: A 26 local C r closing lemma for 3-dimensional systems, Ergodic Theory and Dynamical Systems 27 11 (1990), 443–454. 28 [16] J. Hofbauer and J.W.H. So, Multiple limit cycles for three dimensional LotkaVolterra equa- 29 tions, Appl. Math. Lett. 7 (1994), 65-70. 30 [17] Y.A. Kuznetsov, “Elements of Applied Bifurcation Theory”, Springer-Verlag, 3rd edition, 31 2004. 32 [18] G. Laval and R. Pellat, “Plasma Physics. Proceedings of Summer School of Theoretical 33 Physics”, Gordon and Breach, NY, 1975. 34 [19] J, Llibre, D.D. Novaes and M.A. Teixeira, Higher order averaging theory for finding periodic 35 solutions via Brouwer degree, Nonlinearity 27 (2014), 563–583. 36 [20] J. Llibre and D. Xiao, Limit cycles bifurcating from a non-isolated zero-Hopf equilibrium of 37 3-dimensional differential systems, Proc. Amer. Math. Soc. 142 (2014), 2047–2062. 38 [21] N.G. Lloyd, Degree Theory, Cambridge University Press, 1978. 39 [22] A.J. Lotka, Analytical note on certain rhythmic relations in organic systems, Proc. Natl. 40 Acad. Sci. U.S. 6 (1920), 410–415. 41 [23] Z. Lu and Y. Luo, Two limit cycles in three-dimensional LotkaVolterra systems, Comput. 42 Math. Appl. 44 (2002), 51-66. 43 [24] Z. Lu and Y. Luo, Three limit cycles for a three-dimensional LotkaVolterra competitive 44 system with a heteroclinic cycle, Comput. Math. Appl. 46 (2003), 231-238. 45 [25] R.M. May, “Stability and Complexity in Model Ecosystems”, Princeton, NJ, 1974. 46 [26] J. Scheurle and J. Marsden, Bifurcation to quasi-periodic tori in the interaction of steady 47 state and Hopf bifurcations, SIAM. J. Math. Anal. 15 (1984), 1055–1074. 48 [27] S. Smale, On the differential equations of species in competition, J. Math. Biology 3 (1976), 49 5–7. 50 [28] V. Volterra, “Lecons sur la Théorie Mathématique de la Lutte pour la vie”, Gauthier Villars, 51 Paris, 1931. 52 [29] P. Yu, M. Han and D. Xiao, Four small limit cycles around a Hopf singular point in 3- 53 dimensional competitive LotkaVolterra systems, J. Math. Anal. Appl. 436 (2016), 521–555. 54 [30] D. Xiao and W. Li, Limit cycles for the competitive three dimensional LotkaVolterra system, 55 J. Differential Equations 164 (2000), 1-15. 56 [31] M.L. Zeeman, Hopf bifurcations in competitive three-dimensional Lotka-Volterra systems, 57 Dyn. Stab. Sys. 8 (1993), 189–217. 14 1 E-mail address: mahan@shnu.edu.cn 2 E-mail address: jllibre@mat.uab.cat 3 E-mail address: ytian22@shnu.edu.cn